Background

Overview

We perform multinomial logistic regressions using the multinom function from the nnet package to predict latent class membership (most likely trajectory assignment) in our best GMM model for each symptom measure (anxiety, GAD-7; depression, PHQ-9; anhedonia, MASQ CD anhedonia subscale).

Description of best model for each symptom measure

See GMM notebooks and pre-registration for a description of how best models were selected. An overview of the model for each symptom measure and a brief explanation for its selection is described here.

Depression (PHQ-9)

A 4 class model was selected.

reasons here. Have inserted a place holder of the 5 class model as this is looking best so far. Replace this image with the figure for the best model when that has been identified

best fit model for trajectories of symptoms of depression over one year of the pandemic

We interpret each class as follows:

Class 1 (green line):

Class 2 (green line):

Class 3 (green line):

Class 4 (green line):

Class 5 (green line):

Anxiety (GAD-7)

A 4 class model was selected.

best fit model for trajectories of symptoms of anxiety over one year of the pandemic

We interpret each class as follows:

Class 1 (green line):

Class 2 (green line):

Class 3 (green line):

Class 4 (green line):

Class 5 (green line):

Anhedonia (MASQ)

A x class model was selected.

reasons here. Have inserted a place holder of the 5 class model as this is looking best so far. Replace this image with the figure for the best model when that has been identified

best fit model for trajectories of anhedonia over one year of the pandemic

We interpret each class as follows:

Class 1 (green line):

Class 2 (green line):

Class 3 (green line):

Class 4 (green line):

Class 5 (green line):

Description of selected predictors

The predictors we include which were significant predictors of change in symptoms of anxiety and/or depression in our sample from before the pandemic to the first measurement pint during the pandemic are:

  • Gender (being female predicted worsening of symtoms)
  • Age (being younger; i.e. <35, predicted worsening of symptoms)
  • Employment status pre pandemic (being unemployed predicted worsening of symptoms)
  • Pre-existing mental health diagnoses (a range of mental health diagnoses predicted worsening of symptoms)

We include a few additional measures which were not significant predictors in our first paper, but which we feel are important to consider here.

  • ethnicity (described in the pre-registration; see note below for more detail)
  • key worker status. This was a predictor in our first paper.although we did not plan to include it in the pre-registration, given the piecewise nature of our best fitting latent growth curve based on national pandemic related events, we consider this to be a potential important explanatory variable for differing responses to lockdowns and restrictions easing.
  • employment change during early pandemic (increased employment or became employed, decreased employment or became unemployed, furloughed, no change). Whilst we did not use this variable in our earlier paper, given the role for employment n change from pre-pandemic to early pandemic, it is reasonable to infer that any change in employment status during the pandemic may be an important predictor of symptom trajectories during the pandemic.

A note on how we will handle ethnicity as a predictor.

We will use the ethnicity variable in two different ways in order to ensure we are making the best possible use of the data we have available for this important predictor. First, we will use each self reported ethnicity category as an individual predictor. Although each individual group may be too small to be well powered to detect small effects, we feel it is important to recognise that experiences of different minoritised ethnic communities are not homogeneous, and these different experiences may have an important impact on the overall trajectory of mental healthy symptoms during the pandemic.By collapsing across different minoritsed ethnic community groups, we may miss important differences in potential risk factors and outcomes associated with any one specific ethnicity or background experience.

We will then perform a second analysis, where we group all minoritised ethnic communities together, comparing them to the visible majority white ethnicity group. Whist this grouping does not respect the different experiences that different communities may experience, it does consider the impact of being a “visible minority” and the potential risk for discrimination that this may bring. This analysis will be better powered to detect smaller effects.

We hope that by performing this analysis in both ways, we will be able to fully consider the possible impact of experiences of all ethnic communities on anxiety, depression and anhedonia trajectories throughout the pandemic.

SOME EDITS NEEDED WHEN FINAL CLASS DECIDED

Anywhere that the code is dependent on the number of classes in dataset

Analyses

set up

read in all functions from function library

read in class probability and predictor data

The format of the data output saved from Cprobabilities in mplus is:

first section = that persons raw data for each time point (as many columns as there are time points, so GAD and PHQ = 17, MASQ = 15)
Second section = Individuals Intercept, Lin1 slope, lin2 slope lin 3 slope, lin4 slope (always 5 columns) Third section = Class intercept, ine1 slope, lin2 slope, lin3 slope, lin4 slope (always 5 columns) Fourth section = class posterior probability estimate for each estimated class (as many columns as there are classes)
Second last column: most likely class membership Last column: ID

read in the predictor data

##Data preparation ### name and select the class and ID variables for multinomial regression use our custom function rename_class_data

merge predictors with mplus class data

Left join using MPlus hashed ID so we only retain those included in the MPlus model stage.

This doesn’t change the numbers (reassuring)

drop people without data on any predictor

this represents people who were not missing on all variables as per MPLus output file, and our true analytic sample.

phq.AllData <- phq.AllData[!apply(is.na(phq.AllData[,4:17]), 1, all),]

gad.AllData <- gad.AllData[!apply(is.na(gad.AllData[,4:17]), 1, all),]

masq.AllData <- masq.AllData[!apply(is.na(masq.AllData[,4:17]), 1, all),]

Descriptives of the relationship between class and predictors

PHQ

tbl_cross(phq.AllData, 
          row=MostLikelyClass, col=age,
          percent = "row")
Characteristic age Total
16-18 19-25 26-35 36-45 46-55 56-65 66-70 71-75 76+ Unknown
MostLikelyClass
1 201 (7.2%) 517 (19%) 585 (21%) 485 (17%) 625 (22%) 293 (11%) 53 (1.9%) 18 (0.6%) 5 (0.2%) 0 (0%) 2,782 (100%)
2 92 (0.4%) 844 (4.1%) 2,389 (12%) 2,772 (14%) 4,250 (21%) 5,333 (26%) 2,483 (12%) 1,728 (8.4%) 561 (2.7%) 7 (<0.1%) 20,459 (100%)
3 125 (2.4%) 632 (12%) 1,225 (24%) 999 (19%) 1,093 (21%) 780 (15%) 179 (3.5%) 87 (1.7%) 38 (0.7%) 0 (0%) 5,158 (100%)
4 148 (4.6%) 535 (17%) 832 (26%) 576 (18%) 598 (19%) 403 (13%) 77 (2.4%) 37 (1.1%) 12 (0.4%) 1 (<0.1%) 3,219 (100%)
Total 566 (1.8%) 2,528 (8.0%) 5,031 (16%) 4,832 (15%) 6,566 (21%) 6,809 (22%) 2,792 (8.8%) 1,870 (5.9%) 616 (1.9%) 8 (<0.1%) 31,618 (100%)
tbl_cross(phq.AllData, 
          row=MostLikelyClass, col=gender,
          percent = "row")
Characteristic gender Total
Male Female Non-binary/Prefer to self define Unknown
MostLikelyClass
1 511 (18%) 2,161 (78%) 84 (3.0%) 26 (0.9%) 2,782 (100%)
2 7,607 (37%) 12,735 (62%) 87 (0.4%) 30 (0.1%) 20,459 (100%)
3 1,020 (20%) 4,001 (78%) 110 (2.1%) 27 (0.5%) 5,158 (100%)
4 510 (16%) 2,619 (81%) 76 (2.4%) 14 (0.4%) 3,219 (100%)
Total 9,648 (31%) 21,516 (68%) 357 (1.1%) 97 (0.3%) 31,618 (100%)
tbl_cross(phq.AllData, 
          row=MostLikelyClass, col=ethnicity,
          percent = "row")
Characteristic ethnicity Total
White Mixed or multiple ethnicity Asian Black Arab Other Unknown
MostLikelyClass
1 2,584 (93%) 86 (3.1%) 30 (1.1%) 9 (0.3%) 1 (<0.1%) 24 (0.9%) 48 (1.7%) 2,782 (100%)
2 17,080 (83%) 224 (1.1%) 173 (0.8%) 55 (0.3%) 3 (<0.1%) 49 (0.2%) 2,875 (14%) 20,459 (100%)
3 4,700 (91%) 112 (2.2%) 51 (1.0%) 17 (0.3%) 1 (<0.1%) 32 (0.6%) 245 (4.7%) 5,158 (100%)
4 2,956 (92%) 66 (2.1%) 34 (1.1%) 9 (0.3%) 2 (<0.1%) 31 (1.0%) 121 (3.8%) 3,219 (100%)
Total 27,320 (86%) 488 (1.5%) 288 (0.9%) 90 (0.3%) 7 (<0.1%) 136 (0.4%) 3,289 (10%) 31,618 (100%)
tbl_cross(phq.AllData, 
          row=MostLikelyClass, col=ethnicity_clustered,
          percent = "row")
Characteristic ethnicity_clustered Total
White Minoritised ethnic community Unknown
MostLikelyClass
1 2,584 (93%) 150 (5.4%) 48 (1.7%) 2,782 (100%)
2 17,080 (83%) 504 (2.5%) 2,875 (14%) 20,459 (100%)
3 4,700 (91%) 213 (4.1%) 245 (4.7%) 5,158 (100%)
4 2,956 (92%) 142 (4.4%) 121 (3.8%) 3,219 (100%)
Total 27,320 (86%) 1,009 (3.2%) 3,289 (10%) 31,618 (100%)
tbl_cross(phq.AllData, 
          row=MostLikelyClass, col=employment_prior_covid,
          percent = "row")
Characteristic employment_prior_covid Total
Employed Retired Student Unemployed Unknown
MostLikelyClass
1 1,275 (46%) 196 (7.0%) 226 (8.1%) 463 (17%) 622 (22%) 2,782 (100%)
2 12,730 (62%) 6,127 (30%) 313 (1.5%) 465 (2.3%) 824 (4.0%) 20,459 (100%)
3 3,311 (64%) 533 (10%) 237 (4.6%) 399 (7.7%) 678 (13%) 5,158 (100%)
4 2,039 (63%) 242 (7.5%) 218 (6.8%) 278 (8.6%) 442 (14%) 3,219 (100%)
Total 19,355 (61%) 7,098 (22%) 994 (3.1%) 1,605 (5.1%) 2,566 (8.1%) 31,618 (100%)
tbl_cross(phq.AllData, 
          row=MostLikelyClass, col=employment_change,
          percent = "row")
Characteristic employment_change Total
Decreased employment Employment not changed Furloughed Increased employment Unknown
MostLikelyClass
1 283 (10%) 1,228 (44%) 251 (9.0%) 152 (5.5%) 868 (31%) 2,782 (100%)
2 1,642 (8.0%) 8,566 (42%) 1,875 (9.2%) 571 (2.8%) 7,805 (38%) 20,459 (100%)
3 522 (10%) 2,332 (45%) 560 (11%) 254 (4.9%) 1,490 (29%) 5,158 (100%)
4 383 (12%) 1,371 (43%) 363 (11%) 142 (4.4%) 960 (30%) 3,219 (100%)
Total 2,830 (9.0%) 13,497 (43%) 3,049 (9.6%) 1,119 (3.5%) 11,123 (35%) 31,618 (100%)
tbl_cross(phq.AllData, 
          row=MostLikelyClass, col=key_worker,
          percent = "row")
Characteristic key_worker Total
Government defined key worker Not a key worker Unknown
MostLikelyClass
1 787 (28%) 1,706 (61%) 289 (10%) 2,782 (100%)
2 7,410 (36%) 12,667 (62%) 382 (1.9%) 20,459 (100%)
3 2,055 (40%) 2,776 (54%) 327 (6.3%) 5,158 (100%)
4 1,236 (38%) 1,720 (53%) 263 (8.2%) 3,219 (100%)
Total 11,488 (36%) 18,869 (60%) 1,261 (4.0%) 31,618 (100%)
tbl_cross(phq.AllData, 
          row=MostLikelyClass, col=anxiety_disorders,
          percent = "row")
Characteristic anxiety_disorders Total
Anxiety disorder No anxiety disorder Unknown
MostLikelyClass
1 2,332 (84%) 366 (13%) 84 (3.0%) 2,782 (100%)
2 6,018 (29%) 13,873 (68%) 568 (2.8%) 20,459 (100%)
3 3,632 (70%) 1,305 (25%) 221 (4.3%) 5,158 (100%)
4 2,317 (72%) 738 (23%) 164 (5.1%) 3,219 (100%)
Total 14,299 (45%) 16,282 (51%) 1,037 (3.3%) 31,618 (100%)
tbl_cross(phq.AllData, 
          row=MostLikelyClass, col=depressive_disorders,
          percent = "row")
Characteristic depressive_disorders Total
Depressive disorder No depressive disorder Unknown
MostLikelyClass
1 2,549 (92%) 227 (8.2%) 6 (0.2%) 2,782 (100%)
2 7,376 (36%) 13,014 (64%) 69 (0.3%) 20,459 (100%)
3 4,181 (81%) 968 (19%) 9 (0.2%) 5,158 (100%)
4 2,689 (84%) 526 (16%) 4 (0.1%) 3,219 (100%)
Total 16,795 (53%) 14,735 (47%) 88 (0.3%) 31,618 (100%)
tbl_cross(phq.AllData, 
          row=MostLikelyClass, col=eating_disorders,
          percent = "row")
Characteristic eating_disorders Total
Eating disorder No eating disorder Unknown
MostLikelyClass
1 678 (24%) 2,095 (75%) 9 (0.3%) 2,782 (100%)
2 705 (3.4%) 19,742 (96%) 12 (<0.1%) 20,459 (100%)
3 686 (13%) 4,458 (86%) 14 (0.3%) 5,158 (100%)
4 486 (15%) 2,730 (85%) 3 (<0.1%) 3,219 (100%)
Total 2,555 (8.1%) 29,025 (92%) 38 (0.1%) 31,618 (100%)
tbl_cross(phq.AllData, 
          row=MostLikelyClass, col=obsessive_compulsive_disorders,
          percent = "row")
Characteristic obsessive_compulsive_disorders Total
No obsessive compulsive disorder Obsessive compulsive disorder Unknown
MostLikelyClass
1 2,148 (77%) 632 (23%) 2 (<0.1%) 2,782 (100%)
2 19,797 (97%) 606 (3.0%) 56 (0.3%) 20,459 (100%)
3 4,536 (88%) 618 (12%) 4 (<0.1%) 5,158 (100%)
4 2,778 (86%) 437 (14%) 4 (0.1%) 3,219 (100%)
Total 29,259 (93%) 2,293 (7.3%) 66 (0.2%) 31,618 (100%)
tbl_cross(phq.AllData, 
          row=MostLikelyClass, col=psychotic_disorders,
          percent = "row")
Characteristic psychotic_disorders Total
No psychotic disorder Psychotic disorder Unknown
MostLikelyClass
1 2,618 (94%) 153 (5.5%) 11 (0.4%) 2,782 (100%)
2 20,215 (99%) 173 (0.8%) 71 (0.3%) 20,459 (100%)
3 5,002 (97%) 133 (2.6%) 23 (0.4%) 5,158 (100%)
4 3,144 (98%) 67 (2.1%) 8 (0.2%) 3,219 (100%)
Total 30,979 (98%) 526 (1.7%) 113 (0.4%) 31,618 (100%)
tbl_cross(phq.AllData, 
          row=MostLikelyClass, col=personality_disorders,
          percent = "row")
Characteristic personality_disorders Total
No personality disorder Personality disorder Unknown
MostLikelyClass
1 2,275 (82%) 496 (18%) 11 (0.4%) 2,782 (100%)
2 20,172 (99%) 216 (1.1%) 71 (0.3%) 20,459 (100%)
3 4,783 (93%) 352 (6.8%) 23 (0.4%) 5,158 (100%)
4 2,956 (92%) 255 (7.9%) 8 (0.2%) 3,219 (100%)
Total 30,186 (95%) 1,319 (4.2%) 113 (0.4%) 31,618 (100%)
tbl_cross(phq.AllData, 
          row=MostLikelyClass, col=bipolar_disorders,
          percent = "row")
Characteristic bipolar_disorders Total
Bipolar disorder No bipolar disorder Unknown
MostLikelyClass
1 275 (9.9%) 2,505 (90%) 2 (<0.1%) 2,782 (100%)
2 330 (1.6%) 20,073 (98%) 56 (0.3%) 20,459 (100%)
3 316 (6.1%) 4,838 (94%) 4 (<0.1%) 5,158 (100%)
4 188 (5.8%) 3,027 (94%) 4 (0.1%) 3,219 (100%)
Total 1,109 (3.5%) 30,443 (96%) 66 (0.2%) 31,618 (100%)

GAD

tbl_cross(gad.AllData, 
          row=MostLikelyClass, col=age,
          percent = "row")
Characteristic age Total
16-18 19-25 26-35 36-45 46-55 56-65 66-70 71-75 76+ Unknown
MostLikelyClass
1 190 (0.8%) 1,111 (4.9%) 2,935 (13%) 3,138 (14%) 4,722 (21%) 5,696 (25%) 2,558 (11%) 1,759 (7.8%) 574 (2.5%) 6 (<0.1%) 22,689 (100%)
2 221 (5.7%) 742 (19%) 936 (24%) 694 (18%) 783 (20%) 388 (10%) 67 (1.7%) 24 (0.6%) 10 (0.3%) 2 (<0.1%) 3,867 (100%)
3 99 (4.3%) 356 (15%) 590 (26%) 435 (19%) 426 (18%) 297 (13%) 63 (2.7%) 28 (1.2%) 13 (0.6%) 0 (0%) 2,307 (100%)
4 43 (2.0%) 245 (11%) 464 (21%) 459 (21%) 495 (23%) 350 (16%) 84 (3.8%) 42 (1.9%) 16 (0.7%) 0 (0%) 2,198 (100%)
Total 553 (1.8%) 2,454 (7.9%) 4,925 (16%) 4,726 (15%) 6,426 (21%) 6,731 (22%) 2,772 (8.9%) 1,853 (6.0%) 613 (2.0%) 8 (<0.1%) 31,061 (100%)
tbl_cross(gad.AllData, 
          row=MostLikelyClass, col=gender,
          percent = "row")
Characteristic gender Total
Male Female Non-binary/Prefer to self define Unknown
MostLikelyClass
1 8,111 (36%) 14,365 (63%) 162 (0.7%) 51 (0.2%) 22,689 (100%)
2 592 (15%) 3,132 (81%) 116 (3.0%) 27 (0.7%) 3,867 (100%)
3 347 (15%) 1,917 (83%) 37 (1.6%) 6 (0.3%) 2,307 (100%)
4 407 (19%) 1,740 (79%) 40 (1.8%) 11 (0.5%) 2,198 (100%)
Total 9,457 (30%) 21,154 (68%) 355 (1.1%) 95 (0.3%) 31,061 (100%)
tbl_cross(gad.AllData, 
          row=MostLikelyClass, col=ethnicity,
          percent = "row")
Characteristic ethnicity Total
White Mixed or multiple ethnicity Asian Black Arab Other Unknown
MostLikelyClass
1 19,204 (85%) 285 (1.3%) 193 (0.9%) 59 (0.3%) 4 (<0.1%) 64 (0.3%) 2,880 (13%) 22,689 (100%)
2 3,602 (93%) 96 (2.5%) 38 (1.0%) 9 (0.2%) 1 (<0.1%) 34 (0.9%) 87 (2.2%) 3,867 (100%)
3 2,075 (90%) 54 (2.3%) 25 (1.1%) 8 (0.3%) 2 (<0.1%) 27 (1.2%) 116 (5.0%) 2,307 (100%)
4 2,001 (91%) 43 (2.0%) 25 (1.1%) 14 (0.6%) 0 (0%) 11 (0.5%) 104 (4.7%) 2,198 (100%)
Total 26,882 (87%) 478 (1.5%) 281 (0.9%) 90 (0.3%) 7 (<0.1%) 136 (0.4%) 3,187 (10%) 31,061 (100%)
tbl_cross(gad.AllData, 
          row=MostLikelyClass, col=ethnicity_clustered,
          percent = "row")
Characteristic ethnicity_clustered Total
White Minoritised ethnic community Unknown
MostLikelyClass
1 19,204 (85%) 605 (2.7%) 2,880 (13%) 22,689 (100%)
2 3,602 (93%) 178 (4.6%) 87 (2.2%) 3,867 (100%)
3 2,075 (90%) 116 (5.0%) 116 (5.0%) 2,307 (100%)
4 2,001 (91%) 93 (4.2%) 104 (4.7%) 2,198 (100%)
Total 26,882 (87%) 992 (3.2%) 3,187 (10%) 31,061 (100%)
tbl_cross(gad.AllData, 
          row=MostLikelyClass, col=employment_prior_covid,
          percent = "row")
Characteristic employment_prior_covid Total
Employed Retired Student Unemployed Unknown
MostLikelyClass
1 13,986 (62%) 6,349 (28%) 436 (1.9%) 758 (3.3%) 1,160 (5.1%) 22,689 (100%)
2 2,075 (54%) 250 (6.5%) 298 (7.7%) 531 (14%) 713 (18%) 3,867 (100%)
3 1,510 (65%) 185 (8.0%) 145 (6.3%) 152 (6.6%) 315 (14%) 2,307 (100%)
4 1,474 (67%) 247 (11%) 107 (4.9%) 147 (6.7%) 223 (10%) 2,198 (100%)
Total 19,045 (61%) 7,031 (23%) 986 (3.2%) 1,588 (5.1%) 2,411 (7.8%) 31,061 (100%)
tbl_cross(gad.AllData, 
          row=MostLikelyClass, col=employment_change,
          percent = "row")
Characteristic employment_change Total
Decreased employment Employment not changed Furloughed Increased employment Unknown
MostLikelyClass
1 1,886 (8.3%) 9,545 (42%) 2,100 (9.3%) 667 (2.9%) 8,491 (37%) 22,689 (100%)
2 419 (11%) 1,671 (43%) 383 (9.9%) 230 (5.9%) 1,164 (30%) 3,867 (100%)
3 264 (11%) 1,027 (45%) 258 (11%) 107 (4.6%) 651 (28%) 2,307 (100%)
4 230 (10%) 1,043 (47%) 254 (12%) 94 (4.3%) 577 (26%) 2,198 (100%)
Total 2,799 (9.0%) 13,286 (43%) 2,995 (9.6%) 1,098 (3.5%) 10,883 (35%) 31,061 (100%)
tbl_cross(gad.AllData, 
          row=MostLikelyClass, col=key_worker,
          percent = "row")
Characteristic key_worker Total
Government defined key worker Not a key worker Unknown
MostLikelyClass
1 8,184 (36%) 13,965 (62%) 540 (2.4%) 22,689 (100%)
2 1,271 (33%) 2,245 (58%) 351 (9.1%) 3,867 (100%)
3 925 (40%) 1,217 (53%) 165 (7.2%) 2,307 (100%)
4 930 (42%) 1,198 (55%) 70 (3.2%) 2,198 (100%)
Total 11,310 (36%) 18,625 (60%) 1,126 (3.6%) 31,061 (100%)
tbl_cross(gad.AllData, 
          row=MostLikelyClass, col=anxiety_disorders,
          percent = "row")
Characteristic anxiety_disorders Total
Anxiety disorder No anxiety disorder Unknown
MostLikelyClass
1 7,435 (33%) 14,491 (64%) 763 (3.4%) 22,689 (100%)
2 3,312 (86%) 474 (12%) 81 (2.1%) 3,867 (100%)
3 1,694 (73%) 511 (22%) 102 (4.4%) 2,307 (100%)
4 1,635 (74%) 483 (22%) 80 (3.6%) 2,198 (100%)
Total 14,076 (45%) 15,959 (51%) 1,026 (3.3%) 31,061 (100%)
tbl_cross(gad.AllData, 
          row=MostLikelyClass, col=depressive_disorders,
          percent = "row")
Characteristic depressive_disorders Total
Depressive disorder No depressive disorder Unknown
MostLikelyClass
1 9,595 (42%) 13,028 (57%) 66 (0.3%) 22,689 (100%)
2 3,360 (87%) 499 (13%) 8 (0.2%) 3,867 (100%)
3 1,841 (80%) 463 (20%) 3 (0.1%) 2,307 (100%)
4 1,764 (80%) 433 (20%) 1 (<0.1%) 2,198 (100%)
Total 16,560 (53%) 14,423 (46%) 78 (0.3%) 31,061 (100%)
tbl_cross(gad.AllData, 
          row=MostLikelyClass, col=eating_disorders,
          percent = "row")
Characteristic eating_disorders Total
Eating disorder No eating disorder Unknown
MostLikelyClass
1 1,024 (4.5%) 21,644 (95%) 21 (<0.1%) 22,689 (100%)
2 858 (22%) 3,002 (78%) 7 (0.2%) 3,867 (100%)
3 350 (15%) 1,954 (85%) 3 (0.1%) 2,307 (100%)
4 281 (13%) 1,913 (87%) 4 (0.2%) 2,198 (100%)
Total 2,513 (8.1%) 28,513 (92%) 35 (0.1%) 31,061 (100%)
tbl_cross(gad.AllData, 
          row=MostLikelyClass, col=obsessive_compulsive_disorders,
          percent = "row")
Characteristic obsessive_compulsive_disorders Total
No obsessive compulsive disorder Obsessive compulsive disorder Unknown
MostLikelyClass
1 21,781 (96%) 855 (3.8%) 53 (0.2%) 22,689 (100%)
2 3,049 (79%) 814 (21%) 4 (0.1%) 3,867 (100%)
3 1,991 (86%) 314 (14%) 2 (<0.1%) 2,307 (100%)
4 1,922 (87%) 276 (13%) 0 (0%) 2,198 (100%)
Total 28,743 (93%) 2,259 (7.3%) 59 (0.2%) 31,061 (100%)
tbl_cross(gad.AllData, 
          row=MostLikelyClass, col=psychotic_disorders,
          percent = "row")
Characteristic psychotic_disorders Total
No psychotic disorder Psychotic disorder Unknown
MostLikelyClass
1 22,383 (99%) 234 (1.0%) 72 (0.3%) 22,689 (100%)
2 3,692 (95%) 163 (4.2%) 12 (0.3%) 3,867 (100%)
3 2,248 (97%) 51 (2.2%) 8 (0.3%) 2,307 (100%)
4 2,120 (96%) 70 (3.2%) 8 (0.4%) 2,198 (100%)
Total 30,443 (98%) 518 (1.7%) 100 (0.3%) 31,061 (100%)
tbl_cross(gad.AllData, 
          row=MostLikelyClass, col=personality_disorders,
          percent = "row")
Characteristic personality_disorders Total
No personality disorder Personality disorder Unknown
MostLikelyClass
1 22,203 (98%) 414 (1.8%) 72 (0.3%) 22,689 (100%)
2 3,302 (85%) 553 (14%) 12 (0.3%) 3,867 (100%)
3 2,138 (93%) 161 (7.0%) 8 (0.3%) 2,307 (100%)
4 2,024 (92%) 166 (7.6%) 8 (0.4%) 2,198 (100%)
Total 29,667 (96%) 1,294 (4.2%) 100 (0.3%) 31,061 (100%)
tbl_cross(gad.AllData, 
          row=MostLikelyClass, col=bipolar_disorders,
          percent = "row")
Characteristic bipolar_disorders Total
Bipolar disorder No bipolar disorder Unknown
MostLikelyClass
1 513 (2.3%) 22,123 (98%) 53 (0.2%) 22,689 (100%)
2 321 (8.3%) 3,542 (92%) 4 (0.1%) 3,867 (100%)
3 123 (5.3%) 2,182 (95%) 2 (<0.1%) 2,307 (100%)
4 135 (6.1%) 2,063 (94%) 0 (0%) 2,198 (100%)
Total 1,092 (3.5%) 29,910 (96%) 59 (0.2%) 31,061 (100%)

MASQ

tbl_cross(masq.AllData, 
          row=MostLikelyClass, col=age,
          percent = "row")
Characteristic age Total
16-18 19-25 26-35 36-45 46-55 56-65 66-70 71-75 76+ Unknown
MostLikelyClass
1 8 (1.7%) 32 (6.9%) 89 (19%) 91 (20%) 100 (22%) 90 (19%) 28 (6.0%) 18 (3.9%) 6 (1.3%) 1 (0.2%) 463 (100%)
2 45 (0.6%) 333 (4.7%) 835 (12%) 859 (12%) 1,319 (18%) 1,861 (26%) 988 (14%) 693 (9.7%) 214 (3.0%) 2 (<0.1%) 7,149 (100%)
3 289 (1.7%) 1,187 (7.0%) 2,604 (15%) 2,772 (16%) 3,806 (22%) 3,788 (22%) 1,384 (8.1%) 907 (5.3%) 297 (1.7%) 3 (<0.1%) 17,037 (100%)
4 7 (1.1%) 55 (8.8%) 101 (16%) 95 (15%) 138 (22%) 137 (22%) 61 (9.8%) 22 (3.5%) 8 (1.3%) 0 (0%) 624 (100%)
5 12 (2.4%) 43 (8.4%) 96 (19%) 72 (14%) 109 (21%) 109 (21%) 38 (7.5%) 21 (4.1%) 9 (1.8%) 0 (0%) 509 (100%)
Total 361 (1.4%) 1,650 (6.4%) 3,725 (14%) 3,889 (15%) 5,472 (21%) 5,985 (23%) 2,499 (9.7%) 1,661 (6.4%) 534 (2.1%) 6 (<0.1%) 25,782 (100%)
tbl_cross(masq.AllData, 
          row=MostLikelyClass, col=gender,
          percent = "row")
Characteristic gender Total
Male Female Non-binary/Prefer to self define Unknown
MostLikelyClass
1 127 (27%) 329 (71%) 6 (1.3%) 1 (0.2%) 463 (100%)
2 3,128 (44%) 3,966 (55%) 45 (0.6%) 10 (0.1%) 7,149 (100%)
3 4,481 (26%) 12,278 (72%) 220 (1.3%) 58 (0.3%) 17,037 (100%)
4 120 (19%) 499 (80%) 2 (0.3%) 3 (0.5%) 624 (100%)
5 99 (19%) 400 (79%) 7 (1.4%) 3 (0.6%) 509 (100%)
Total 7,955 (31%) 17,472 (68%) 280 (1.1%) 75 (0.3%) 25,782 (100%)
tbl_cross(masq.AllData, 
          row=MostLikelyClass, col=ethnicity,
          percent = "row")
Characteristic ethnicity Total
White Mixed or multiple ethnicity Asian Black Arab Other Unknown
MostLikelyClass
1 403 (87%) 9 (1.9%) 4 (0.9%) 1 (0.2%) 0 (0%) 3 (0.6%) 43 (9.3%) 463 (100%)
2 5,892 (82%) 82 (1.1%) 65 (0.9%) 22 (0.3%) 1 (<0.1%) 15 (0.2%) 1,072 (15%) 7,149 (100%)
3 15,152 (89%) 254 (1.5%) 124 (0.7%) 48 (0.3%) 5 (<0.1%) 87 (0.5%) 1,367 (8.0%) 17,037 (100%)
4 524 (84%) 11 (1.8%) 9 (1.4%) 2 (0.3%) 0 (0%) 1 (0.2%) 77 (12%) 624 (100%)
5 451 (89%) 10 (2.0%) 7 (1.4%) 0 (0%) 0 (0%) 5 (1.0%) 36 (7.1%) 509 (100%)
Total 22,422 (87%) 366 (1.4%) 209 (0.8%) 73 (0.3%) 6 (<0.1%) 111 (0.4%) 2,595 (10%) 25,782 (100%)
tbl_cross(masq.AllData, 
          row=MostLikelyClass, col=ethnicity_clustered,
          percent = "row")
Characteristic ethnicity_clustered Total
White Minoritised ethnic community Unknown
MostLikelyClass
1 403 (87%) 17 (3.7%) 43 (9.3%) 463 (100%)
2 5,892 (82%) 185 (2.6%) 1,072 (15%) 7,149 (100%)
3 15,152 (89%) 518 (3.0%) 1,367 (8.0%) 17,037 (100%)
4 524 (84%) 23 (3.7%) 77 (12%) 624 (100%)
5 451 (89%) 22 (4.3%) 36 (7.1%) 509 (100%)
Total 22,422 (87%) 765 (3.0%) 2,595 (10%) 25,782 (100%)
tbl_cross(masq.AllData, 
          row=MostLikelyClass, col=employment_prior_covid,
          percent = "row")
Characteristic employment_prior_covid Total
Employed Retired Student Unemployed Unknown
MostLikelyClass
1 331 (71%) 75 (16%) 16 (3.5%) 23 (5.0%) 18 (3.9%) 463 (100%)
2 4,386 (61%) 2,335 (33%) 110 (1.5%) 139 (1.9%) 179 (2.5%) 7,149 (100%)
3 10,544 (62%) 3,694 (22%) 582 (3.4%) 1,109 (6.5%) 1,108 (6.5%) 17,037 (100%)
4 417 (67%) 141 (23%) 28 (4.5%) 14 (2.2%) 24 (3.8%) 624 (100%)
5 351 (69%) 103 (20%) 17 (3.3%) 25 (4.9%) 13 (2.6%) 509 (100%)
Total 16,029 (62%) 6,348 (25%) 753 (2.9%) 1,310 (5.1%) 1,342 (5.2%) 25,782 (100%)
tbl_cross(masq.AllData, 
          row=MostLikelyClass, col=employment_change,
          percent = "row")
Characteristic employment_change Total
Decreased employment Employment not changed Furloughed Increased employment Unknown
MostLikelyClass
1 54 (12%) 225 (49%) 44 (9.5%) 14 (3.0%) 126 (27%) 463 (100%)
2 554 (7.7%) 2,915 (41%) 668 (9.3%) 201 (2.8%) 2,811 (39%) 7,149 (100%)
3 1,637 (9.6%) 7,659 (45%) 1,601 (9.4%) 644 (3.8%) 5,496 (32%) 17,037 (100%)
4 68 (11%) 248 (40%) 71 (11%) 23 (3.7%) 214 (34%) 624 (100%)
5 69 (14%) 199 (39%) 63 (12%) 23 (4.5%) 155 (30%) 509 (100%)
Total 2,382 (9.2%) 11,246 (44%) 2,447 (9.5%) 905 (3.5%) 8,802 (34%) 25,782 (100%)
tbl_cross(masq.AllData, 
          row=MostLikelyClass, col=key_worker,
          percent = "row")
Characteristic key_worker Total
Government defined key worker Not a key worker Unknown
MostLikelyClass
1 186 (40%) 277 (60%) 0 (0%) 463 (100%)
2 2,569 (36%) 4,536 (63%) 44 (0.6%) 7,149 (100%)
3 6,295 (37%) 10,551 (62%) 191 (1.1%) 17,037 (100%)
4 253 (41%) 366 (59%) 5 (0.8%) 624 (100%)
5 209 (41%) 299 (59%) 1 (0.2%) 509 (100%)
Total 9,512 (37%) 16,029 (62%) 241 (0.9%) 25,782 (100%)
tbl_cross(masq.AllData, 
          row=MostLikelyClass, col=anxiety_disorders,
          percent = "row")
Characteristic anxiety_disorders Total
Anxiety disorder No anxiety disorder Unknown
MostLikelyClass
1 243 (52%) 192 (41%) 28 (6.0%) 463 (100%)
2 1,716 (24%) 5,290 (74%) 143 (2.0%) 7,149 (100%)
3 8,875 (52%) 7,529 (44%) 633 (3.7%) 17,037 (100%)
4 259 (42%) 342 (55%) 23 (3.7%) 624 (100%)
5 253 (50%) 232 (46%) 24 (4.7%) 509 (100%)
Total 11,346 (44%) 13,585 (53%) 851 (3.3%) 25,782 (100%)
tbl_cross(masq.AllData, 
          row=MostLikelyClass, col=depressive_disorders,
          percent = "row")
Characteristic depressive_disorders Total
Depressive disorder No depressive disorder Unknown
MostLikelyClass
1 305 (66%) 157 (34%) 1 (0.2%) 463 (100%)
2 2,002 (28%) 5,124 (72%) 23 (0.3%) 7,149 (100%)
3 10,571 (62%) 6,436 (38%) 30 (0.2%) 17,037 (100%)
4 303 (49%) 320 (51%) 1 (0.2%) 624 (100%)
5 318 (62%) 190 (37%) 1 (0.2%) 509 (100%)
Total 13,499 (52%) 12,227 (47%) 56 (0.2%) 25,782 (100%)
tbl_cross(masq.AllData, 
          row=MostLikelyClass, col=eating_disorders,
          percent = "row")
Characteristic eating_disorders Total
Eating disorder No eating disorder Unknown
MostLikelyClass
1 42 (9.1%) 420 (91%) 1 (0.2%) 463 (100%)
2 222 (3.1%) 6,922 (97%) 5 (<0.1%) 7,149 (100%)
3 1,569 (9.2%) 15,453 (91%) 15 (<0.1%) 17,037 (100%)
4 46 (7.4%) 577 (92%) 1 (0.2%) 624 (100%)
5 52 (10%) 457 (90%) 0 (0%) 509 (100%)
Total 1,931 (7.5%) 23,829 (92%) 22 (<0.1%) 25,782 (100%)
tbl_cross(masq.AllData, 
          row=MostLikelyClass, col=obsessive_compulsive_disorders,
          percent = "row")
Characteristic obsessive_compulsive_disorders Total
No obsessive compulsive disorder Obsessive compulsive disorder Unknown
MostLikelyClass
1 440 (95%) 22 (4.8%) 1 (0.2%) 463 (100%)
2 6,914 (97%) 213 (3.0%) 22 (0.3%) 7,149 (100%)
3 15,564 (91%) 1,454 (8.5%) 19 (0.1%) 17,037 (100%)
4 586 (94%) 37 (5.9%) 1 (0.2%) 624 (100%)
5 483 (95%) 25 (4.9%) 1 (0.2%) 509 (100%)
Total 23,987 (93%) 1,751 (6.8%) 44 (0.2%) 25,782 (100%)
tbl_cross(masq.AllData, 
          row=MostLikelyClass, col=psychotic_disorders,
          percent = "row")
Characteristic psychotic_disorders Total
No psychotic disorder Psychotic disorder Unknown
MostLikelyClass
1 457 (99%) 5 (1.1%) 1 (0.2%) 463 (100%)
2 7,064 (99%) 60 (0.8%) 25 (0.3%) 7,149 (100%)
3 16,662 (98%) 330 (1.9%) 45 (0.3%) 17,037 (100%)
4 608 (97%) 12 (1.9%) 4 (0.6%) 624 (100%)
5 496 (97%) 13 (2.6%) 0 (0%) 509 (100%)
Total 25,287 (98%) 420 (1.6%) 75 (0.3%) 25,782 (100%)
tbl_cross(masq.AllData, 
          row=MostLikelyClass, col=personality_disorders,
          percent = "row")
Characteristic personality_disorders Total
No personality disorder Personality disorder Unknown
MostLikelyClass
1 439 (95%) 23 (5.0%) 1 (0.2%) 463 (100%)
2 7,061 (99%) 63 (0.9%) 25 (0.3%) 7,149 (100%)
3 16,136 (95%) 856 (5.0%) 45 (0.3%) 17,037 (100%)
4 604 (97%) 16 (2.6%) 4 (0.6%) 624 (100%)
5 486 (95%) 23 (4.5%) 0 (0%) 509 (100%)
Total 24,726 (96%) 981 (3.8%) 75 (0.3%) 25,782 (100%)
tbl_cross(masq.AllData, 
          row=MostLikelyClass, col=bipolar_disorders,
          percent = "row")
Characteristic bipolar_disorders Total
Bipolar disorder No bipolar disorder Unknown
MostLikelyClass
1 28 (6.0%) 434 (94%) 1 (0.2%) 463 (100%)
2 91 (1.3%) 7,036 (98%) 22 (0.3%) 7,149 (100%)
3 712 (4.2%) 16,306 (96%) 19 (0.1%) 17,037 (100%)
4 16 (2.6%) 607 (97%) 1 (0.2%) 624 (100%)
5 25 (4.9%) 483 (95%) 1 (0.2%) 509 (100%)
Total 872 (3.4%) 24,866 (96%) 44 (0.2%) 25,782 (100%)

specify reference category for predictors

make age midpoint the reference for age, no change in employment the reference for employment change, all disorders have the no disorder group as reference, not a keyworker reference for key worker status

# PHQ
phq.AllData$age <- factor(phq.AllData$age,
                                       levels=c("16-18","19-25","26-35","36-45","46-55","56-65","66-70","71-75","76+")) %>%
  relevel(phq.AllData$age, ref = "26-35")

phq.AllData$employment_change <- factor(phq.AllData$employment_change,
                                       levels=c("Decreased employment","Employment not changed","Furloughed","Increased employment")) %>%
  relevel(phq.AllData$employment_change, ref = "Employment not changed")

phq.AllData$key_worker <- factor(phq.AllData$key_worker,
                                       levels=c("Government defined key worker","Not a key worker")) %>%
  relevel(phq.AllData$key_worker, ref = "Not a key worker")

phq.AllData$anxiety_disorders <- factor(phq.AllData$anxiety_disorders,
                                       levels=c("Anxiety disorder","No anxiety disorder")) %>%
  relevel(phq.AllData$anxiety_disorders, ref = "No anxiety disorder")

phq.AllData$depressive_disorders <- factor(phq.AllData$depressive_disorders,
                                       levels=c("Depressive disorder","No depressive disorder")) %>%
  relevel(phq.AllData$depressive_disorders, ref = "No depressive disorder")

phq.AllData$eating_disorders <- factor(phq.AllData$eating_disorders,
                                       levels=c("Eating disorder","No eating disorder")) %>%
  relevel(phq.AllData$eating_disorders, ref = "No eating disorder")

phq.AllData$obsessive_compulsive_disorders <- factor(phq.AllData$obsessive_compulsive_disorders,
                                       levels=c("Obsessive compulsive disorder","No obsessive compulsive disorder")) %>%
  relevel(phq.AllData$obsessive_compulsive_disorders, ref = "No obsessive compulsive disorder")

phq.AllData$psychotic_disorders <- factor(phq.AllData$psychotic_disorders,
                                       levels=c("Psychotic disorder","No psychotic disorder")) %>%
  relevel(phq.AllData$psychotic_disorders, ref = "No psychotic disorder")


phq.AllData$personality_disorders <- factor(phq.AllData$personality_disorders,
                                       levels=c("Personality disorder","No personality disorder")) %>%
  relevel(phq.AllData$personality_disorders, ref = "No personality disorder")


phq.AllData$bipolar_disorders <- factor(phq.AllData$bipolar_disorders,
                                       levels=c("Bipolar disorder","No bipolar disorder")) %>%
  relevel(phq.AllData$bipolar_disorders, ref = "No bipolar disorder")




# GAD
gad.AllData$age <- factor(gad.AllData$age,
                                       levels=c("16-18","19-25","26-35","36-45","46-55","56-65","66-70","71-75","76+")) %>%
  relevel(gad.AllData$age, ref = "26-35")

gad.AllData$employment_change <- factor(gad.AllData$employment_change,
                                       levels=c("Decreased employment","Employment not changed","Furloughed","Increased employment")) %>%
  relevel(gad.AllData$employment_change, ref = "Employment not changed")

gad.AllData$key_worker <- factor(gad.AllData$key_worker,
                                       levels=c("Government defined key worker","Not a key worker")) %>%
  relevel(gad.AllData$key_worker, ref = "Not a key worker")

gad.AllData$anxiety_disorders <- factor(gad.AllData$anxiety_disorders,
                                       levels=c("Anxiety disorder","No anxiety disorder")) %>%
  relevel(gad.AllData$anxiety_disorders, ref = "No anxiety disorder")

gad.AllData$depressive_disorders <- factor(gad.AllData$depressive_disorders,
                                       levels=c("Depressive disorder","No depressive disorder")) %>%
  relevel(gad.AllData$depressive_disorders, ref = "No depressive disorder")

gad.AllData$eating_disorders <- factor(gad.AllData$eating_disorders,
                                       levels=c("Eating disorder","No eating disorder")) %>%
  relevel(gad.AllData$eating_disorders, ref = "No eating disorder")

gad.AllData$obsessive_compulsive_disorders <- factor(gad.AllData$obsessive_compulsive_disorders,
                                       levels=c("Obsessive compulsive disorder","No obsessive compulsive disorder")) %>%
  relevel(gad.AllData$obsessive_compulsive_disorders, ref = "No obsessive compulsive disorder")

gad.AllData$psychotic_disorders <- factor(gad.AllData$psychotic_disorders,
                                       levels=c("Psychotic disorder","No psychotic disorder")) %>%
  relevel(gad.AllData$psychotic_disorders, ref = "No psychotic disorder")


gad.AllData$personality_disorders <- factor(gad.AllData$personality_disorders,
                                       levels=c("Personality disorder","No personality disorder")) %>%
  relevel(gad.AllData$personality_disorders, ref = "No personality disorder")


gad.AllData$bipolar_disorders <- factor(gad.AllData$bipolar_disorders,
                                       levels=c("Bipolar disorder","No bipolar disorder")) %>%
  relevel(gad.AllData$bipolar_disorders, ref = "No bipolar disorder")



# MASQ
masq.AllData$age <- factor(masq.AllData$age,
                                       levels=c("16-18","19-25","26-35","36-45","46-55","56-65","66-70","71-75","76+")) %>%
  relevel(masq.AllData$age, ref = "26-35")

masq.AllData$employment_change <- factor(masq.AllData$employment_change,
                                       levels=c("Decreased employment","Employment not changed","Furloughed","Increased employment")) %>%
  relevel(masq.AllData$employment_change, ref = "Employment not changed")

masq.AllData$key_worker <- factor(masq.AllData$key_worker,
                                       levels=c("Government defined key worker","Not a key worker")) %>%
  relevel(masq.AllData$key_worker, ref = "Not a key worker")

masq.AllData$anxiety_disorders <- factor(masq.AllData$anxiety_disorders,
                                       levels=c("Anxiety disorder","No anxiety disorder")) %>%
  relevel(masq.AllData$anxiety_disorders, ref = "No anxiety disorder")

masq.AllData$depressive_disorders <- factor(masq.AllData$depressive_disorders,
                                       levels=c("Depressive disorder","No depressive disorder")) %>%
  relevel(masq.AllData$depressive_disorders, ref = "No depressive disorder")

masq.AllData$eating_disorders <- factor(masq.AllData$eating_disorders,
                                       levels=c("Eating disorder","No eating disorder")) %>%
  relevel(masq.AllData$eating_disorders, ref = "No eating disorder")

masq.AllData$obsessive_compulsive_disorders <- factor(masq.AllData$obsessive_compulsive_disorders,
                                       levels=c("Obsessive compulsive disorder","No obsessive compulsive disorder")) %>%
  relevel(masq.AllData$obsessive_compulsive_disorders, ref = "No obsessive compulsive disorder")

masq.AllData$psychotic_disorders <- factor(masq.AllData$psychotic_disorders,
                                       levels=c("Psychotic disorder","No psychotic disorder")) %>%
  relevel(masq.AllData$psychotic_disorders, ref = "No psychotic disorder")


masq.AllData$personality_disorders <- factor(masq.AllData$personality_disorders,
                                       levels=c("Personality disorder","No personality disorder")) %>%
  relevel(masq.AllData$personality_disorders, ref = "No personality disorder")


masq.AllData$bipolar_disorders <- factor(masq.AllData$bipolar_disorders,
                                       levels=c("Bipolar disorder","No bipolar disorder")) %>%
  relevel(masq.AllData$bipolar_disorders, ref = "No bipolar disorder")

Multinomial regression

formulas

the formulas will be the same for all disorders so create these once here

formula_clustered <- "MostLikelyClass ~ age + gender + ethnicity_clustered + employment_prior_covid + employment_change +key_worker + anxiety_disorders + depressive_disorders + eating_disorders +obsessive_compulsive_disorders +psychotic_disorders + personality_disorders + bipolar_disorders"

formula <- "MostLikelyClass ~ age + gender + ethnicity + employment_prior_covid + employment_change +key_worker + anxiety_disorders + depressive_disorders + eating_disorders +obsessive_compulsive_disorders +psychotic_disorders + personality_disorders + bipolar_disorders"

PHQ

prepare and run models

specify reference category for outcome

will use the largest class group as reference

BiggestClass.phq <- phq.AllData %>%
  group_by(MostLikelyClass) %>%
  summarise(n=n()) %>%
  arrange(-n) %>%
  slice(1)%>%
  select(MostLikelyClass) %>%
  as.character(.)

BiggestClass.phq
[1] "2"

Will need to adapt this to the correct number of classes depending on final model

phq.AllData$MostLikelyClass <- factor(phq.AllData$MostLikelyClass,
                                       levels=c("1","2","3","4")) %>%
  relevel(phq.AllData$MostLikelyClass, ref = BiggestClass.phq)
Perform regression (ethnicity clustered)
phq.model_clustered <- multinom(formula_clustered,
                            data = phq.AllData)
# weights:  108 (78 variable)
initial  value 22776.816353 
iter  10 value 15413.454463
iter  20 value 15018.010656
iter  30 value 14829.289484
iter  40 value 14788.839015
iter  50 value 14780.869242
iter  60 value 14778.570012
iter  70 value 14777.735248
iter  80 value 14777.582866
iter  90 value 14777.495686
final  value 14777.495046 
converged
Perform regression (ethnicity as is)
phq.model<- multinom(formula,
                            data = phq.AllData)
# weights:  124 (90 variable)
initial  value 22776.816353 
iter  10 value 15411.839789
iter  20 value 15013.647665
iter  30 value 14829.869090
iter  40 value 14785.994207
iter  50 value 14778.358654
iter  60 value 14775.657187
iter  70 value 14774.939734
iter  80 value 14774.798621
iter  90 value 14774.630785
iter 100 value 14774.597139
final  value 14774.597139 
stopped after 100 iterations

View and process results

1: coefficients 2: std errors 3: residual deviances & AIC ##### view summary: Ethnicity clustered

summary(phq.model_clustered)
Call:
multinom(formula = formula_clustered, data = phq.AllData)

Coefficients:
  (Intercept)  age16-18  age19-25   age36-45   age46-55   age56-65  age66-70
1   -4.549853 2.1048558 0.7790551 -0.3149913 -0.2498514 -0.9543816 -2.341322
3   -2.547683 0.9017894 0.3711199 -0.2213567 -0.3897026 -0.7070398 -1.121361
4   -3.140364 1.3448804 0.5455621 -0.3948331 -0.6598507 -1.0175320 -1.234575
   age71-75     age76+ genderFemale genderNon-binary/Prefer to self define
1 -1.019484 -10.730111  -0.02542089                              0.5807275
3 -1.453851  -1.269494   0.16522410                              0.5840882
4 -2.661366  -1.382587   0.30206855                              0.7226346
  ethnicity_clusteredMinoritised ethnic community employment_prior_covidRetired
1                                       0.3708296                     1.7704825
3                                       0.2395683                     0.6350815
4                                       0.2208408                     1.2337246
  employment_prior_covidStudent employment_prior_covidUnemployed
1                     1.8104805                         1.960040
3                     0.4844898                         1.036247
4                     0.8510254                         1.242283
  employment_changeDecreased employment employment_changeFurloughed
1                             0.5650230                   0.4106201
3                             0.3115720                   0.2948073
4                             0.5195315                   0.2970980
  employment_changeIncreased employment key_workerGovernment defined key worker
1                            0.13575139                               0.2522181
3                            0.19257965                               0.1636542
4                            0.02743296                               0.1640294
  anxiety_disordersAnxiety disorder depressive_disordersDepressive disorder
1                         0.9085664                                1.948273
3                         0.6100999                                1.209666
4                         0.5927057                                1.287837
  eating_disordersEating disorder
1                       0.9211555
3                       0.5823767
4                       0.5287970
  obsessive_compulsive_disordersObsessive compulsive disorder
1                                                   0.5989446
3                                                   0.2161699
4                                                   0.3800703
  psychotic_disordersPsychotic disorder
1                            -0.1139183
3                            -0.1722263
4                            -0.5338968
  personality_disordersPersonality disorder bipolar_disordersBipolar disorder
1                                 1.0722802                         0.6964569
3                                 0.4741771                         0.5781367
4                                 0.4548947                         0.5946864

Std. Errors:
  (Intercept)  age16-18   age19-25   age36-45   age46-55   age56-65  age66-70
1  0.14816774 0.2316826 0.10736351 0.09418929 0.08985454 0.12326540 0.6038183
3  0.08507648 0.2312577 0.08919550 0.06657359 0.06493759 0.07797673 0.2208516
4  0.10673324 0.2304602 0.09590673 0.07867265 0.07868899 0.09853262 0.2723463
   age71-75       age76+ genderFemale genderNon-binary/Prefer to self define
1 0.4966426 7.368934e-06   0.08179817                              0.2451186
3 0.3554065 6.205729e-01   0.05707060                              0.2163430
4 0.7222626 7.493670e-01   0.07267529                              0.2381598
  ethnicity_clusteredMinoritised ethnic community employment_prior_covidRetired
1                                       0.1488697                     0.3823456
3                                       0.1139642                     0.3232712
4                                       0.1338969                     0.3511517
  employment_prior_covidStudent employment_prior_covidUnemployed
1                     0.4933996                       0.10826932
3                     0.5260395                       0.09334759
4                     0.5513469                       0.10525114
  employment_changeDecreased employment employment_changeFurloughed
1                            0.09267171                  0.09352993
3                            0.06770352                  0.06549133
4                            0.07816144                  0.07892212
  employment_changeIncreased employment key_workerGovernment defined key worker
1                            0.12813110                              0.07569931
3                            0.09580487                              0.05127901
4                            0.11850412                              0.06199946
  anxiety_disordersAnxiety disorder depressive_disordersDepressive disorder
1                        0.09314783                              0.11596604
3                        0.05791695                              0.06135151
4                        0.07176675                              0.07758567
  eating_disordersEating disorder
1                      0.09503508
3                      0.08237463
4                      0.09344419
  obsessive_compulsive_disordersObsessive compulsive disorder
1                                                  0.09393152
3                                                  0.08383057
4                                                  0.09228084
  psychotic_disordersPsychotic disorder
1                             0.1921388
3                             0.1726906
4                             0.2178653
  personality_disordersPersonality disorder bipolar_disordersBipolar disorder
1                                 0.1274057                         0.1391339
3                                 0.1236615                         0.1230321
4                                 0.1371931                         0.1400993

Residual Deviance: 29554.99 
AIC: 29710.99 
view summary: Ethnicity as is
summary(phq.model)
Call:
multinom(formula = formula, data = phq.AllData)

Coefficients:
  (Intercept)  age16-18  age19-25   age36-45   age46-55   age56-65  age66-70
1   -4.553307 2.1039961 0.7796502 -0.3135582 -0.2467409 -0.9515374 -2.336530
3   -2.547118 0.9002638 0.3705796 -0.2201190 -0.3894447 -0.7063024 -1.119318
4   -3.134719 1.3428491 0.5459123 -0.3963309 -0.6621591 -1.0204405 -1.239464
   age71-75     age76+ genderFemale genderNon-binary/Prefer to self define
1 -1.017143 -12.224352  -0.02480251                              0.5807562
3 -1.453370  -1.270669   0.16369397                              0.5807002
4 -2.662583  -1.382245   0.30174999                              0.7205345
  ethnicityMixed or multiple ethnicity ethnicityAsian ethnicityBlack
1                            0.4629838     0.40858586      0.2066424
3                            0.3457871     0.12366733      0.1582022
4                            0.1690055     0.08441411      0.3413721
  ethnicityArab ethnicityOther employment_prior_covidRetired
1     -5.662052      0.1570573                     1.7674622
3      1.029187      0.1020597                     0.6327436
4     -7.525484      0.5425931                     1.2358961
  employment_prior_covidStudent employment_prior_covidUnemployed
1                     1.8080856                         1.961047
3                     0.4822341                         1.037672
4                     0.8539063                         1.240262
  employment_changeDecreased employment employment_changeFurloughed
1                             0.5664661                   0.4110791
3                             0.3117479                   0.2944403
4                             0.5180671                   0.2954592
  employment_changeIncreased employment key_workerGovernment defined key worker
1                            0.13930698                               0.2508914
3                            0.19402322                               0.1642079
4                            0.02349589                               0.1640721
  anxiety_disordersAnxiety disorder depressive_disordersDepressive disorder
1                         0.9106050                                1.947505
3                         0.6097613                                1.209428
4                         0.5902519                                1.286131
  eating_disordersEating disorder
1                       0.9218347
3                       0.5834878
4                       0.5270773
  obsessive_compulsive_disordersObsessive compulsive disorder
1                                                   0.5979472
3                                                   0.2154201
4                                                   0.3810171
  psychotic_disordersPsychotic disorder
1                            -0.1152364
3                            -0.1729216
4                            -0.5311309
  personality_disordersPersonality disorder bipolar_disordersBipolar disorder
1                                 1.0726504                         0.6981503
3                                 0.4744518                         0.5786918
4                                 0.4551671                         0.5916305

Std. Errors:
  (Intercept)  age16-18   age19-25   age36-45   age46-55   age56-65  age66-70
1  0.14831580 0.2316582 0.10738478 0.09423458 0.08991647 0.12331345 0.6036621
3  0.08517708 0.2312299 0.08922588 0.06658966 0.06497050 0.07801082 0.2208580
4  0.10682035 0.2305134 0.09593152 0.07869585 0.07872043 0.09856665 0.2724340
   age71-75       age76+ genderFemale genderNon-binary/Prefer to self define
1 0.4967691 1.749907e-06   0.08183394                              0.2452735
3 0.3554055 6.209694e-01   0.05709798                              0.2163849
4 0.7215908 7.480872e-01   0.07270660                              0.2382940
  ethnicityMixed or multiple ethnicity ethnicityAsian ethnicityBlack
1                            0.2016425      0.2955720      0.5253943
3                            0.1572987      0.2192610      0.3860519
4                            0.1909652      0.2666643      0.4347993
  ethnicityArab ethnicityOther employment_prior_covidRetired
1     49.416999      0.3884027                     0.3823841
3      1.416853      0.3144694                     0.3233035
4     85.478333      0.3207641                     0.3511165
  employment_prior_covidStudent employment_prior_covidUnemployed
1                     0.4934754                        0.1082854
3                     0.5261723                        0.0933672
4                     0.5513879                        0.1053026
  employment_changeDecreased employment employment_changeFurloughed
1                            0.09268309                  0.09356790
3                            0.06770987                  0.06551691
4                            0.07818836                  0.07894781
  employment_changeIncreased employment key_workerGovernment defined key worker
1                            0.12822434                              0.07572031
3                            0.09582231                              0.05129569
4                            0.11857777                              0.06202039
  anxiety_disordersAnxiety disorder depressive_disordersDepressive disorder
1                        0.09319813                              0.11601099
3                        0.05796675                              0.06137905
4                        0.07182591                              0.07762486
  eating_disordersEating disorder
1                      0.09507111
3                      0.08239786
4                      0.09345317
  obsessive_compulsive_disordersObsessive compulsive disorder
1                                                  0.09395176
3                                                  0.08385171
4                                                  0.09229920
  psychotic_disordersPsychotic disorder
1                             0.1921904
3                             0.1727261
4                             0.2178664
  personality_disordersPersonality disorder bipolar_disordersBipolar disorder
1                                 0.1274059                         0.1391611
3                                 0.1236762                         0.1230710
4                                 0.1372178                         0.1401488

Residual Deviance: 29549.19 
AIC: 29729.19 
create z and p values

make Z

z_clustered.phq <- summary(phq.model_clustered)$coefficients/summary(phq.model_clustered)$standard.errors
  
z_as.is.phq <- summary(phq.model)$coefficients/summary(phq.model)$standard.errors

print z clusted ethnicity model

z_clustered.phq
  (Intercept) age16-18 age19-25  age36-45  age46-55   age56-65  age66-70
1   -30.70745 9.085084 7.256237 -3.344237 -2.780620  -7.742493 -3.877527
3   -29.94580 3.899499 4.160747 -3.324993 -6.001187  -9.067317 -5.077441
4   -29.42255 5.835631 5.688466 -5.018682 -8.385554 -10.326855 -4.533106
   age71-75        age76+ genderFemale genderNon-binary/Prefer to self define
1 -2.052753 -1.456128e+06   -0.3107757                               2.369170
3 -4.090670 -2.045681e+00    2.8950826                               2.699825
4 -3.684763 -1.845007e+00    4.1564134                               3.034242
  ethnicity_clusteredMinoritised ethnic community employment_prior_covidRetired
1                                        2.490967                      4.630581
3                                        2.102137                      1.964547
4                                        1.649335                      3.513366
  employment_prior_covidStudent employment_prior_covidUnemployed
1                     3.6693996                         18.10338
3                     0.9210142                         11.10095
4                     1.5435391                         11.80303
  employment_changeDecreased employment employment_changeFurloughed
1                              6.097038                    4.390253
3                              4.602005                    4.501470
4                              6.646903                    3.764445
  employment_changeIncreased employment key_workerGovernment defined key worker
1                             1.0594727                                3.331842
3                             2.0101239                                3.191447
4                             0.2314938                                2.645659
  anxiety_disordersAnxiety disorder depressive_disordersDepressive disorder
1                          9.754026                                16.80038
3                         10.534048                                19.71697
4                          8.258779                                16.59891
  eating_disordersEating disorder
1                        9.692795
3                        7.069855
4                        5.658961
  obsessive_compulsive_disordersObsessive compulsive disorder
1                                                    6.376396
3                                                    2.578652
4                                                    4.118626
  psychotic_disordersPsychotic disorder
1                            -0.5928962
3                            -0.9973114
4                            -2.4505819
  personality_disordersPersonality disorder bipolar_disordersBipolar disorder
1                                  8.416263                          5.005661
3                                  3.834477                          4.699070
4                                  3.315726                          4.244750

print z ethnicity as is

z_as.is.phq
  (Intercept) age16-18 age19-25  age36-45  age46-55   age56-65  age66-70
1   -30.70008 9.082330 7.260341 -3.327422 -2.744113  -7.716412 -3.870594
3   -29.90379 3.893371 4.153275 -3.305603 -5.994177  -9.053903 -5.068045
4   -29.34571 5.825472 5.690645 -5.036236 -8.411528 -10.352797 -4.549595
   age71-75        age76+ genderFemale genderNon-binary/Prefer to self define
1 -2.047517 -6.985717e+06   -0.3030834                               2.367790
3 -4.089330 -2.046267e+00    2.8668961                               2.683645
4 -3.689879 -1.847706e+00    4.1502423                               3.023720
  ethnicityMixed or multiple ethnicity ethnicityAsian ethnicityBlack
1                            2.2960627      1.3823566      0.3933091
3                            2.1982830      0.5640188      0.4097952
4                            0.8850066      0.3165557      0.7851257
  ethnicityArab ethnicityOther employment_prior_covidRetired
1   -0.11457701      0.4043671                      4.622216
3    0.72638972      0.3245459                      1.957120
4   -0.08803967      1.6915639                      3.519903
  employment_prior_covidStudent employment_prior_covidUnemployed
1                     3.6639837                         18.10999
3                     0.9164946                         11.11388
4                     1.5486490                         11.77808
  employment_changeDecreased employment employment_changeFurloughed
1                              6.111861                    4.393378
3                              4.604172                    4.494111
4                              6.625886                    3.742463
  employment_changeIncreased employment key_workerGovernment defined key worker
1                             1.0864317                                3.313397
3                             2.0248230                                3.201204
4                             0.1981475                                2.645454
  anxiety_disordersAnxiety disorder depressive_disordersDepressive disorder
1                          9.770635                                16.78725
3                         10.519155                                19.70425
4                          8.217813                                16.56855
  eating_disordersEating disorder
1                        9.696265
3                        7.081347
4                        5.640016
  obsessive_compulsive_disordersObsessive compulsive disorder
1                                                    6.364407
3                                                    2.569060
4                                                    4.128065
  psychotic_disordersPsychotic disorder
1                            -0.5995952
3                            -1.0011317
4                            -2.4378744
  personality_disordersPersonality disorder bipolar_disordersBipolar disorder
1                                  8.419158                          5.016848
3                                  3.836242                          4.702097
4                                  3.317113                          4.221445

two tailed z test p values

p_clustered.phq <- (1 - pnorm(abs(z_clustered.phq), 0, 1)) * 2
p_as.is.phq <- (1 - pnorm(abs(z_as.is.phq), 0, 1)) * 2

show p for two tailed z test for ethnicity clustered model

p_clustered.phq 
  (Intercept)     age16-18     age19-25     age36-45     age46-55     age56-65
1           0 0.000000e+00 3.979039e-13 8.250917e-04 5.425514e-03 9.769963e-15
3           0 9.639175e-05 3.172090e-05 8.842073e-04 1.958799e-09 0.000000e+00
4           0 5.358735e-09 1.281860e-08 5.202710e-07 0.000000e+00 0.000000e+00
      age66-70     age71-75     age76+ genderFemale
1 1.055237e-04 4.009657e-02 0.00000000 7.559711e-01
3 3.825528e-07 4.301282e-05 0.04078774 3.790588e-03
4 5.812267e-06 2.289157e-04 0.06503655 3.232826e-05
  genderNon-binary/Prefer to self define
1                            0.017828077
3                            0.006937603
4                            0.002411408
  ethnicity_clusteredMinoritised ethnic community employment_prior_covidRetired
1                                      0.01273958                  3.646404e-06
3                                      0.03554124                  4.946666e-02
4                                      0.09907909                  4.424671e-04
  employment_prior_covidStudent employment_prior_covidUnemployed
1                  0.0002431208                                0
3                  0.3570430153                                0
4                  0.1227000329                                0
  employment_changeDecreased employment employment_changeFurloughed
1                          1.080516e-09                1.132186e-05
3                          4.184431e-06                6.748506e-06
4                          2.993250e-11                1.669193e-04
  employment_changeIncreased employment key_workerGovernment defined key worker
1                            0.28938457                            0.0008627318
3                            0.04441808                            0.0014156211
4                            0.81693124                            0.0081532077
  anxiety_disordersAnxiety disorder depressive_disordersDepressive disorder
1                      0.000000e+00                                       0
3                      0.000000e+00                                       0
4                      2.220446e-16                                       0
  eating_disordersEating disorder
1                    0.000000e+00
3                    1.550982e-12
4                    1.522920e-08
  obsessive_compulsive_disordersObsessive compulsive disorder
1                                                1.813041e-10
3                                                9.918667e-03
4                                                3.811378e-05
  psychotic_disordersPsychotic disorder
1                            0.55325065
3                            0.31861340
4                            0.01426255
  personality_disordersPersonality disorder bipolar_disordersBipolar disorder
1                              0.0000000000                      5.567075e-07
3                              0.0001258319                      2.613485e-06
4                              0.0009140554                      2.188370e-05

show p for two tailed z test for ethnicity as is model

p_as.is.phq 
  (Intercept)     age16-18     age19-25     age36-45     age46-55     age56-65
1           0 0.000000e+00 3.861356e-13 8.765361e-04 6.067473e-03 1.199041e-14
3           0 9.886074e-05 3.277503e-05 9.477232e-04 2.045184e-09 0.000000e+00
4           0 5.695158e-09 1.265603e-08 4.747745e-07 0.000000e+00 0.000000e+00
      age66-70     age71-75     age76+ genderFemale
1 1.085707e-04 4.060738e-02 0.00000000 7.618263e-01
3 4.019231e-07 4.326218e-05 0.04073007 4.145191e-03
4 5.374926e-06 2.243604e-04 0.06464496 3.321235e-05
  genderNon-binary/Prefer to self define ethnicityMixed or multiple ethnicity
1                            0.017894695                           0.02167230
3                            0.007282448                           0.02792895
4                            0.002496871                           0.37615309
  ethnicityAsian ethnicityBlack ethnicityArab ethnicityOther
1      0.1668622      0.6940912     0.9087804     0.68594278
3      0.5727413      0.6819562     0.4675999     0.74552479
4      0.7515808      0.4323799     0.9298451     0.09072915
  employment_prior_covidRetired employment_prior_covidStudent
1                  3.796619e-06                  0.0002483226
3                  5.033338e-02                  0.3594075444
4                  4.317043e-04                  0.1214661249
  employment_prior_covidUnemployed employment_changeDecreased employment
1                                0                          9.847620e-10
3                                0                          4.141100e-06
4                                0                          3.451728e-11
  employment_changeFurloughed employment_changeIncreased employment
1                1.116030e-05                            0.27728806
3                6.986125e-06                            0.04288554
4                1.822257e-04                            0.84292964
  key_workerGovernment defined key worker anxiety_disordersAnxiety disorder
1                            0.0009217019                      0.000000e+00
3                            0.0013685480                      0.000000e+00
4                            0.0081581445                      2.220446e-16
  depressive_disordersDepressive disorder eating_disordersEating disorder
1                                       0                    0.000000e+00
3                                       0                    1.427525e-12
4                                       0                    1.700348e-08
  obsessive_compulsive_disordersObsessive compulsive disorder
1                                                1.960461e-10
3                                                1.019747e-02
4                                                3.658293e-05
  psychotic_disordersPsychotic disorder
1                            0.54877604
3                            0.31676312
4                            0.01477391
  personality_disordersPersonality disorder bipolar_disordersBipolar disorder
1                              0.0000000000                      5.252597e-07
3                              0.0001249311                      2.575029e-06
4                              0.0009095269                      2.427412e-05

Relative risk

Relative risk is the ratio of the probability of choosing one outcome category over the probability of choosing the baseline category.

Relative risk of 1 means there is no difference between outcome and baseline. Greater than 1 indicates outcome is more likely than baseline. Less than one indicates it is less likely than baseline.

RR 1.4 would mean there is a 1.4% chance of outcome occurring relative to baseline etc

if something is below one (e.g. 0.6 x the risk) you can also interpret this as 40% less probability of outcomes (1-0.6 = 0.4). So, if women are 0.6 x as likely to belong to a class as men, that means that there is a 40% lower probability of being in that category for women.

exponentiate the model coefficient to obtain relative risk for all predictors. Can see equations here

clustered ethnicity
relative_risk_clustered.phq <- data.frame(exp(coef(phq.model_clustered))) %>%
  rownames_to_column() %>%
  rename(class = "rowname",
    intercept = "X.Intercept.") %>%
  pivot_longer(!class) %>%
  rename(predictor = name,
         RelativeRisk = value) 


relative_risk_clustered.phq$predictor <- gsub("76.","76+",relative_risk_clustered.phq$predictor)


relative_risk_clustered.phq 
ethnicity
relative_risk_as.is.phq <-data.frame(exp(coef(phq.model))) %>%
  rownames_to_column() %>%
  rename(class = "rowname",
    intercept = "X.Intercept.") %>%
  pivot_longer(!class) %>%
  rename(predictor = name,
         RelativeRisk = value)


relative_risk_as.is.phq$predictor <- gsub("76.","76+",relative_risk_as.is.phq$predictor)

head(relative_risk_as.is.phq)

confidence intervals (rr scale)

The confint package has a specific method that is run for multinom objects (see methods(confint)) Obtain the exponent of upper and lower to get relative risk version

First create the CI (converted to relative risk) and reformat it so it is ready to merge with relative risk data for plotting

clustered ethnicity

CI_clustered.phq <- data.frame(exp(confint(phq.model_clustered,level=0.95))) %>%
  rownames_to_column()

names(CI_clustered.phq) <- c("predictor",
                              "1_lower","1_upper",
                              "3_lower","3_upper",
                              "4_lower","4_upper"
                              )

CI_clustered.phq <- CI_clustered.phq %>%
  pivot_longer(!predictor) %>%
  separate(name, c("class","CI")) %>%
  pivot_wider(names_from = CI,
              values_from= value) 


CI_clustered.phq$predictor <- gsub("-",".",CI_clustered.phq$predictor)
CI_clustered.phq$predictor <- gsub(" ",".",CI_clustered.phq$predictor)
CI_clustered.phq$predictor <- gsub("/",".",CI_clustered.phq$predictor)
CI_clustered.phq$predictor <- gsub("\\(Intercept\\)","intercept",CI_clustered.phq$predictor)


  

head(CI_clustered.phq)

ethnicity as is

CI_as.is.phq<- data.frame(exp(confint(phq.model,level=0.95))) %>%
  rownames_to_column()

names(CI_as.is.phq) <-  c("predictor",
                              "1_lower","1_upper",
                              "3_lower","3_upper",
                              "4_lower","4_upper"
                              )

CI_as.is.phq <- CI_as.is.phq %>%
  pivot_longer(!predictor) %>%
  separate(name, c("class","CI")) %>%
  pivot_wider(names_from = CI,
              values_from= value)

CI_as.is.phq$predictor <- gsub("-",".",CI_as.is.phq$predictor)
CI_as.is.phq$predictor <- gsub(" ",".",CI_as.is.phq$predictor)
CI_as.is.phq$predictor <- gsub("/",".",CI_as.is.phq$predictor)
CI_as.is.phq$predictor <- gsub("\\(Intercept\\)","intercept",CI_as.is.phq$predictor)

head(CI_as.is.phq)

forest plots of relative risk

PHQ

prepare and format plot data

plot_phq_clustered_dat$category <-plot_phq_clustered_dat$predictor

plot_phq_clustered_dat <-
  plot_phq_clustered_dat %>%
  mutate(category = case_when(grepl("age",predictor) ~ "Age",
                              grepl("gender",predictor) ~ "Gender",
                              grepl("ethnicity_clustered",predictor) ~ "Ethnicity",
                              grepl("employment_prior_covid",predictor) ~ "Employment before covid",
                              grepl("employment_change",predictor) ~ "Change in employment due to covid",
                              grepl("key_worker",predictor) ~ "Key worker status",
                              grepl("anxiety_disorders",predictor) ~ "Pre-existing mental health diagnosis",
                              grepl("depressive_disorders",predictor) ~ "Pre-existing mental health diagnosis",
                              grepl("eating_disorders",predictor) ~ "Pre-existing mental health diagnosis",
                              grepl("obsessive_compulsive_disorders",predictor) ~ "Pre-existing mental health diagnosis",
                              grepl("psychotic_disorders",predictor) ~ "Pre-existing mental health diagnosis",
                              grepl("personality_disorders",predictor) ~ "Pre-existing mental health diagnosis",
                              grepl("bipolar_disorders",predictor) ~ "Pre-existing mental health diagnosis",
                              
                              grepl("intercept",predictor) ~ "Intercept"
                              ),
         
         predictor = case_when(grepl("age",predictor) ~ gsub("age","",predictor),
                               grepl("gender",predictor) ~ gsub("gender","",predictor),
                               grepl("ethnicity_clustered",predictor) ~ gsub("ethnicity_clustered","",predictor),
                               grepl("employment_prior_covid",predictor) ~ gsub("employment_prior_covid","",predictor),
                               grepl("employment_change",predictor) ~ gsub("employment_change","",predictor),
                               grepl("key_worker",predictor) ~ gsub("key_worker","",predictor),
                               grepl("anxiety_disorders",predictor) ~ gsub("anxiety_disorders","",predictor),
                               grepl("depressive_disorders",predictor) ~ gsub("depressive_disorders","",predictor),
                               grepl("eating_disorders",predictor) ~ gsub("eating_disorders","",predictor),
                               grepl("obsessive_compulsive_disorders",predictor) ~ gsub("obsessive_compulsive_disorders","",predictor),
                               grepl("psychotic_disorders",predictor) ~ gsub("psychotic_disorders","",predictor),
                               grepl("personality_disorders",predictor) ~ gsub("personality_disorders","",predictor),
                               grepl("bipolar_disorders",predictor) ~ gsub("bipolar_disorders","",predictor),
                               grepl("intercept",predictor) ~ "intercept")
                               ) 
         

plot_phq_as.is_dat$category <-plot_phq_as.is_dat$predictor

plot_phq_as.is_dat <-
plot_phq_as.is_dat %>%
  mutate(category = case_when(grepl("age",predictor) ~ "Age",
                              grepl("gender",predictor) ~ "Gender",
                              grepl("ethnicity",predictor) ~ "Ethnicity",
                              grepl("employment_prior_covid",predictor) ~ "Employment before covid",
                              grepl("employment_change",predictor) ~ "Change in employment due to covid",
                              grepl("key_worker",predictor) ~ "Key worker status",
                              grepl("anxiety_disorders",predictor) ~ "Pre-existing mental health diagnosis",
                              grepl("depressive_disorders",predictor) ~ "Pre-existing mental health diagnosis",
                              grepl("eating_disorders",predictor) ~ "Pre-existing mental health diagnosis",
                              grepl("obsessive_compulsive_disorders",predictor) ~ "Pre-existing mental health diagnosis",
                              grepl("psychotic_disorders",predictor) ~ "Pre-existing mental health diagnosis",
                              grepl("personality_disorders",predictor) ~ "Pre-existing mental health diagnosis",
                              grepl("bipolar_disorders",predictor) ~ "Pre-existing mental health diagnosis",
                              grepl("intercept",predictor) ~ "Intercept"
                              ),
         
         predictor = case_when(grepl("age",predictor) ~ gsub("age","",predictor),
                               grepl("gender",predictor) ~ gsub("gender","",predictor),
                               grepl("ethnicity",predictor) ~ gsub("ethnicity_clustered","",predictor),
                               grepl("employment_prior_covid",predictor) ~ gsub("employment_prior_covid","",predictor),
                               grepl("employment_change",predictor) ~ gsub("employment_change","",predictor),
                               grepl("key_worker",predictor) ~ gsub("key_worker","",predictor),
                               grepl("anxiety_disorders",predictor) ~ gsub("anxiety_disorders","",predictor),
                               grepl("depressive_disorders",predictor) ~ gsub("depressive_disorders","",predictor),
                               grepl("eating_disorders",predictor) ~ gsub("eating_disorders","",predictor),
                               grepl("obsessive_compulsive_disorders",predictor) ~ gsub("obsessive_compulsive_disorders","",predictor),
                               grepl("psychotic_disorders",predictor) ~ gsub("psychotic_disorders","",predictor),
                               grepl("personality_disorders",predictor) ~ gsub("personality_disorders","",predictor),
                               grepl("bipolar_disorders",predictor) ~ gsub("bipolar_disorders","",predictor),
                               grepl("intercept",predictor) ~ "intercept")
                               )

sort variables so categories are together

plot_phq_clustered_dat$category <- factor(plot_phq_clustered_dat$category,
                                       levels=c("Intercept","Gender","Ethnicity","Age","Employment before covid","Change in employment due to covid","Key worker status","Pre-existing mental health diagnosis"))

plot_phq_clustered_dat$predictor <- fct_reorder2(plot_phq_clustered_dat$predictor, plot_phq_clustered_dat$class,  plot_phq_clustered_dat$category)




plot_phq_as.is_dat$category <- factor(plot_phq_as.is_dat$category,
                                       levels=c("Intercept","Gender","Ethnicity","Age","Employment before covid","Change in employment due to covid","Key worker status","Pre-existing mental health diagnosis"))

plot_phq_as.is_dat$predictor <- fct_reorder2(plot_phq_as.is_dat$predictor, plot_phq_as.is_dat$class,  plot_phq_as.is_dat$category)

forest plot clustered ethnicity

Note: smaller category error bars extend beyond a reasonable axis. Wil need fixing for publication.

phq_eth_clustered_forest <-relrisk.forest(plot_phq_clustered_dat,
plot_title="Relative risk (95% CI) of class membership: depression",
sub_title="ethnicity clustered",
lowlim = 0.001,
uplim = 10)


phq_eth_clustered_forest
Warning: Removed 1 rows containing missing values (geom_point).

forest plot ethnicity as is

Note: smaller category error bars extend beyond a reasonable axis. Wil need fixing for publication.

phq_eth_as.is_forest <- relrisk.forest(plot_phq_as.is_dat,
plot_title="Relative risk (95% CI) of class membership: depression",
sub_title="ethnicity as is",
lowlim = 0.001,
uplim = 10)

phq_eth_as.is_forest
Warning: Removed 2 rows containing missing values (geom_point).

#### save plots

ggsave(
  filename=file.path(dirname(dirname(getwd())),"output/ForestPlots/PHQ/depression_RelativeRisk_forestPlot_ethnicityClustered.png"),
  plot=phq_eth_clustered_forest, 
  width = 18, height = 10, dpi = 150, units = "in")
Warning: Removed 1 rows containing missing values (geom_point).
ggsave(
  filename=file.path(dirname(dirname(getwd())),"output/ForestPlots/PHQ/depression_RelativeRisk_forestPlot.png"),
  plot=phq_eth_as.is_forest, 
  width = 18, height = 10, dpi = 150, units = "in")
Warning: Removed 2 rows containing missing values (geom_point).

GAD

prepare and run models

specify reference category for outcome

will use the largest class group as reference

BiggestClass.gad <- gad.AllData %>%
  group_by(MostLikelyClass) %>%
  summarise(n=n()) %>%
  arrange(-n) %>%
  slice(1)%>%
  select(MostLikelyClass) %>%
  as.character(.)

BiggestClass.gad
[1] "1"

Will need to adapt this to the correct number of classes depending on final model

gad.AllData$MostLikelyClass <- factor(gad.AllData$MostLikelyClass,
                                       levels=c("1","2","3","4")) %>%
  relevel(gad.AllData$MostLikelyClass, ref = BiggestClass.gad)
Perform regression (ethnicity clustered)
gad.model_clustered <- multinom(formula_clustered,
                            data = gad.AllData)
# weights:  108 (78 variable)
initial  value 22478.763066 
iter  10 value 14824.087869
iter  20 value 13416.842887
iter  30 value 13142.479298
iter  40 value 13071.880866
iter  50 value 13049.903786
iter  60 value 13034.906346
iter  70 value 13032.108909
iter  80 value 13031.900653
iter  90 value 13031.859302
final  value 13031.857624 
converged
Perform regression (ethnicity as is)
gad.model<- multinom(formula,
                            data = gad.AllData)
# weights:  124 (90 variable)
initial  value 22478.763066 
iter  10 value 14813.923252
iter  20 value 13422.343932
iter  30 value 13137.531913
iter  40 value 13064.638921
iter  50 value 13042.619448
iter  60 value 13027.215094
iter  70 value 13024.194017
iter  80 value 13024.022070
iter  90 value 13023.887259
iter 100 value 13023.867221
final  value 13023.867221 
stopped after 100 iterations

View and process results

1: coefficients 2: std errors 3: residual deviances & AIC ##### view summary: Ethnicity clustered

summary(gad.model_clustered)
Call:
multinom(formula = formula_clustered, data = gad.AllData)

Coefficients:
  (Intercept)  age16-18  age19-25   age36-45     age46-55   age56-65  age66-70
2   -3.650424 1.1671214 0.5167330 -0.3085380 -0.343064521 -0.7999582 -2.365594
3   -3.475616 1.0496975 0.3779647 -0.2286434 -0.534428182 -0.7652060 -1.189222
4   -3.673601 0.1425086 0.3313403  0.1073131 -0.005039378 -0.4684441 -0.753368
    age71-75      age76+ genderFemale genderNon-binary/Prefer to self define
2 -0.9020561  -1.4383121    0.2448057                              0.6393261
3 -1.4549421   0.5644221    0.4638501                              0.2542844
4 -0.3070324 -10.1905261    0.3037216                              0.5180137
  ethnicity_clusteredMinoritised ethnic community employment_prior_covidRetired
2                                       0.0592478                     1.0102689
3                                       0.2659093                    -0.3294263
4                                       0.1400527                     0.4950650
  employment_prior_covidStudent employment_prior_covidUnemployed
2                     1.2912962                        1.1014962
3                     0.5522132                        0.3372711
4                     0.4724846                        0.3397331
  employment_changeDecreased employment employment_changeFurloughed
2                             0.4543803                   0.2636938
3                             0.2660162                   0.1221211
4                             0.1787511                   0.2495741
  employment_changeIncreased employment key_workerGovernment defined key worker
2                            0.30259582                               0.1241316
3                            0.06364632                               0.1524750
4                           -0.01003740                               0.1400552
  anxiety_disordersAnxiety disorder depressive_disordersDepressive disorder
2                          1.341481                               0.8639027
3                          1.017514                               0.5559895
4                          1.014918                               0.6965658
  eating_disordersEating disorder
2                       0.6927949
3                       0.4714443
4                       0.3468174
  obsessive_compulsive_disordersObsessive compulsive disorder
2                                                   0.5245066
3                                                   0.3518988
4                                                   0.3286507
  psychotic_disordersPsychotic disorder
2                           -0.27256668
3                           -0.37148039
4                           -0.01097946
  personality_disordersPersonality disorder bipolar_disordersBipolar disorder
2                                 0.7749959                         0.3766943
3                                 0.3466961                         0.3097247
4                                 0.4947044                         0.3160350

Std. Errors:
  (Intercept)  age16-18   age19-25   age36-45   age46-55   age56-65  age66-70
2   0.1115779 0.1821392 0.08662416 0.07599623 0.07379979 0.09749191 0.5171385
3   0.1209820 0.2100379 0.10223948 0.08583027 0.08811655 0.11020444 0.3503438
4   0.1212064 0.3056786 0.11251066 0.08751113 0.08621793 0.11082011 0.3118558
   age71-75       age76+ genderFemale genderNon-binary/Prefer to self define
2 0.4102723 1.034236e+00   0.07006872                              0.1955054
3 0.5916734 4.990777e-01   0.08497578                              0.2801124
4 0.3794280 5.425413e-06   0.07958294                              0.2547380
  ethnicity_clusteredMinoritised ethnic community employment_prior_covidRetired
2                                       0.1296068                     0.3478351
3                                       0.1405756                     0.6066659
4                                       0.1492563                     0.3994340
  employment_prior_covidStudent employment_prior_covidUnemployed
2                     0.4284839                       0.08869572
3                     0.5929406                       0.11846162
4                     0.5930770                       0.11775712
  employment_changeDecreased employment employment_changeFurloughed
2                            0.07550160                  0.07667607
3                            0.08745712                  0.08780095
4                            0.08980712                  0.08564925
  employment_changeIncreased employment key_workerGovernment defined key worker
2                             0.1025110                              0.06079071
3                             0.1267405                              0.06833989
4                             0.1318353                              0.06847648
  anxiety_disordersAnxiety disorder depressive_disordersDepressive disorder
2                        0.07929849                              0.08052156
3                        0.08484729                              0.08492792
4                        0.08406461                              0.08673466
  eating_disordersEating disorder
2                      0.07895395
3                      0.09658515
4                      0.10152638
  obsessive_compulsive_disordersObsessive compulsive disorder
2                                                  0.07843148
3                                                  0.09753485
4                                                  0.09930082
  psychotic_disordersPsychotic disorder
2                             0.1687666
3                             0.2298316
4                             0.2003119
  personality_disordersPersonality disorder bipolar_disordersBipolar disorder
2                                 0.1052235                         0.1170634
3                                 0.1400174                         0.1485819
4                                 0.1360787                         0.1438603

Residual Deviance: 26063.72 
AIC: 26219.72 
view summary: Ethnicity as is
summary(gad.model)
Call:
multinom(formula = formula, data = gad.AllData)

Coefficients:
  (Intercept)  age16-18  age19-25   age36-45    age46-55   age56-65   age66-70
2   -3.646697 1.1667348 0.5156172 -0.3110604 -0.34601962 -0.8031963 -2.3728387
3   -3.469596 1.0509495 0.3793907 -0.2312247 -0.53736659 -0.7691526 -1.1995273
4   -3.686265 0.1473397 0.3303717  0.1070361 -0.00347322 -0.4660158 -0.7482131
    age71-75     age76+ genderFemale genderNon-binary/Prefer to self define
2 -0.9040785 -1.4422427    0.2447487                              0.6397616
3 -1.4554868  0.5597647    0.4647165                              0.2525852
4 -0.3024041 -8.7921182    0.3078017                              0.5318378
  ethnicityMixed or multiple ethnicity ethnicityAsian ethnicityBlack
2                          -0.05176409  -0.0003202253     -0.1430201
3                           0.12618985   0.1578156311      0.3250710
4                          -0.07803545   0.5743700321      0.6327083
  ethnicityArab ethnicityOther employment_prior_covidRetired
2      1.322949      0.4771204                     1.0146878
3     -6.923904      0.7987081                    -0.3221631
4     -6.750052     -0.5325362                     0.5006125
  employment_prior_covidStudent employment_prior_covidUnemployed
2                     1.2991101                        1.1010880
3                     0.5649541                        0.3347179
4                     0.4802050                        0.3397831
  employment_changeDecreased employment employment_changeFurloughed
2                             0.4524217                   0.2638354
3                             0.2650956                   0.1206205
4                             0.1798413                   0.2530216
  employment_changeIncreased employment key_workerGovernment defined key worker
2                            0.30181772                               0.1250909
3                            0.05921025                               0.1526558
4                           -0.01028472                               0.1390436
  anxiety_disordersAnxiety disorder depressive_disordersDepressive disorder
2                          1.338478                               0.8645838
3                          1.014348                               0.5540216
4                          1.019877                               0.7019997
  eating_disordersEating disorder
2                       0.6921110
3                       0.4679514
4                       0.3450005
  obsessive_compulsive_disordersObsessive compulsive disorder
2                                                   0.5257258
3                                                   0.3536987
4                                                   0.3306319
  psychotic_disordersPsychotic disorder
2                           -0.26935131
3                           -0.36799398
4                           -0.01354605
  personality_disordersPersonality disorder bipolar_disordersBipolar disorder
2                                 0.7743597                         0.3760422
3                                 0.3471852                         0.3058703
4                                 0.4947032                         0.3173593

Std. Errors:
  (Intercept)  age16-18   age19-25   age36-45   age46-55  age56-65  age66-70
2   0.1116937 0.1822423 0.08665573 0.07603363 0.07384173 0.0975301 0.5173850
3   0.1211115 0.2101064 0.10226111 0.08587586 0.08816022 0.1102456 0.3507171
4   0.1214361 0.3057609 0.11255457 0.08755993 0.08626089 0.1108650 0.3119076
   age71-75    age76+ genderFemale genderNon-binary/Prefer to self define
2 0.4101659  1.034099   0.07010760                              0.1955637
3 0.5907168  0.498965   0.08504013                              0.2803780
4 0.3795802 50.638438   0.07965851                              0.2548989
  ethnicityMixed or multiple ethnicity ethnicityAsian ethnicityBlack
2                            0.1783553      0.2787059      0.4671697
3                            0.1969998      0.2953030      0.4636567
4                            0.2179014      0.2547151      0.4131507
  ethnicityArab ethnicityOther employment_prior_covidRetired
2      1.425252      0.3005387                     0.3477641
3     72.668274      0.3211487                     0.6060899
4     69.796153      0.5325911                     0.3993792
  employment_prior_covidStudent employment_prior_covidUnemployed
2                     0.4286759                       0.08874296
3                     0.5920381                       0.11854461
4                     0.5934081                       0.11780540
  employment_changeDecreased employment employment_changeFurloughed
2                            0.07555035                  0.07670187
3                            0.08748728                  0.08783670
4                            0.08984714                  0.08571794
  employment_changeIncreased employment key_workerGovernment defined key worker
2                             0.1025941                              0.06083141
3                             0.1268516                              0.06837122
4                             0.1319701                              0.06849874
  anxiety_disordersAnxiety disorder depressive_disordersDepressive disorder
2                        0.07934682                              0.08058216
3                        0.08491949                              0.08497846
4                        0.08418799                              0.08681954
  eating_disordersEating disorder
2                      0.07895700
3                      0.09660906
4                      0.10161878
  obsessive_compulsive_disordersObsessive compulsive disorder
2                                                  0.07843933
3                                                  0.09755387
4                                                  0.09935671
  psychotic_disordersPsychotic disorder
2                             0.1687378
3                             0.2298358
4                             0.2004391
  personality_disordersPersonality disorder bipolar_disordersBipolar disorder
2                                 0.1052320                         0.1171095
3                                 0.1400425                         0.1486702
4                                 0.1360719                         0.1439566

Residual Deviance: 26047.73 
AIC: 26227.73 
create z and p values

make Z

z_clustered.gad <- summary(gad.model_clustered)$coefficients/summary(gad.model_clustered)$standard.errors
  
z_as.is.gad <- summary(gad.model)$coefficients/summary(gad.model)$standard.errors

print z clusted ethnicity model

z_clustered.gad
  (Intercept) age16-18 age19-25  age36-45   age46-55  age56-65  age66-70
2   -32.71636 6.407855 5.965230 -4.059912 -4.6485841 -8.205380 -4.574390
3   -28.72836 4.997658 3.696856 -2.663902 -6.0650148 -6.943513 -3.394444
4   -30.30865 0.466204 2.944968  1.226280 -0.0584493 -4.227068 -2.415758
    age71-75        age76+ genderFemale genderNon-binary/Prefer to self define
2 -2.1986765 -1.390700e+00     3.493794                              3.2701194
3 -2.4590291  1.130930e+00     5.458616                              0.9077944
4 -0.8091979 -1.878295e+06     3.816416                              2.0335153
  ethnicity_clusteredMinoritised ethnic community employment_prior_covidRetired
2                                       0.4571349                      2.904448
3                                       1.8915760                     -0.543011
4                                       0.9383370                      1.239416
  employment_prior_covidStudent employment_prior_covidUnemployed
2                     3.0136401                        12.418819
3                     0.9313128                         2.847091
4                     0.7966666                         2.885032
  employment_changeDecreased employment employment_changeFurloughed
2                              6.018154                    3.439063
3                              3.041675                    1.390886
4                              1.990389                    2.913908
  employment_changeIncreased employment key_workerGovernment defined key worker
2                            2.95183704                                2.041951
3                            0.50217813                                2.231128
4                           -0.07613593                                2.045303
  anxiety_disordersAnxiety disorder depressive_disordersDepressive disorder
2                          16.91685                               10.728838
3                          11.99229                                6.546605
4                          12.07307                                8.030998
  eating_disordersEating disorder
2                        8.774671
3                        4.881126
4                        3.416032
  obsessive_compulsive_disordersObsessive compulsive disorder
2                                                    6.687450
3                                                    3.607929
4                                                    3.309647
  psychotic_disordersPsychotic disorder
2                           -1.61505149
3                           -1.61631565
4                           -0.05481183
  personality_disordersPersonality disorder bipolar_disordersBipolar disorder
2                                  7.365239                          3.217865
3                                  2.476093                          2.084539
4                                  3.635428                          2.196818

print z ethnicity as is

z_as.is.gad
  (Intercept) age16-18 age19-25  age36-45    age46-55  age56-65  age66-70
2   -32.64908 6.402108 5.950181 -4.091089 -4.68596312 -8.235368 -4.586215
3   -28.64795 5.001987 3.710020 -2.692545 -6.09534096 -6.976716 -3.420213
4   -30.35559 0.481879 2.935213  1.222433 -0.04026414 -4.203454 -2.398830
    age71-75     age76+ genderFemale genderNon-binary/Prefer to self define
2 -2.2041775 -1.3946859     3.491043                              3.2713725
3 -2.4639334  1.1218518     5.464674                              0.9008739
4 -0.7966804 -0.1736254     3.864015                              2.0864655
  ethnicityMixed or multiple ethnicity ethnicityAsian ethnicityBlack
2                           -0.2902302   -0.001148972     -0.3061417
3                            0.6405584    0.534419244      0.7011028
4                           -0.3581228    2.254951039      1.5314226
  ethnicityArab ethnicityOther employment_prior_covidRetired
2    0.92822109       1.587550                     2.9177471
3   -0.09528097       2.487035                    -0.5315435
4   -0.09671094      -0.999897                     1.2534765
  employment_prior_covidStudent employment_prior_covidUnemployed
2                     3.0305182                        12.407610
3                     0.9542529                         2.823561
4                     0.8092322                         2.884274
  employment_changeDecreased employment employment_changeFurloughed
2                              5.988347                    3.439752
3                              3.030104                    1.373236
4                              2.001637                    2.951793
  employment_changeIncreased employment key_workerGovernment defined key worker
2                            2.94186230                                2.056354
3                            0.46676803                                2.232750
4                           -0.07793218                                2.029870
  anxiety_disordersAnxiety disorder depressive_disordersDepressive disorder
2                          16.86871                               10.729220
3                          11.94482                                6.519554
4                          12.11428                                8.085734
  eating_disordersEating disorder
2                        8.765671
3                        4.843763
4                        3.395046
  obsessive_compulsive_disordersObsessive compulsive disorder
2                                                    6.702325
3                                                    3.625676
4                                                    3.327726
  psychotic_disordersPsychotic disorder
2                           -1.59627176
3                           -1.60111702
4                           -0.06758187
  personality_disordersPersonality disorder bipolar_disordersBipolar disorder
2                                  7.358595                          3.211031
3                                  2.479142                          2.057374
4                                  3.635602                          2.204548

two tailed z test p values

p_clustered.gad <- (1 - pnorm(abs(z_clustered.gad), 0, 1)) * 2
p_as.is.gad <- (1 - pnorm(abs(z_as.is.gad), 0, 1)) * 2

show p for two tailed z test for ethnicity clustered model

p_clustered.gad 
  (Intercept)     age16-18     age19-25     age36-45     age46-55     age56-65
2           0 1.475813e-10 2.442897e-09 4.909115e-05 3.342214e-06 2.220446e-16
3           0 5.803079e-07 2.182857e-04 7.724014e-03 1.319415e-09 3.824718e-12
4           0 6.410695e-01 3.229883e-03 2.200935e-01 9.533907e-01 2.367561e-05
      age66-70   age71-75    age76+ genderFemale
2 4.776088e-06 0.02790093 0.1643165 4.762076e-04
3 6.876817e-04 0.01393133 0.2580844 4.798605e-08
4 1.570251e-02 0.41840135 0.0000000 1.354044e-04
  genderNon-binary/Prefer to self define
2                            0.001075021
3                            0.363986836
4                            0.042000490
  ethnicity_clusteredMinoritised ethnic community employment_prior_covidRetired
2                                      0.64757408                   0.003679016
3                                      0.05854749                   0.587122199
4                                      0.34807126                   0.215191371
  employment_prior_covidStudent employment_prior_covidUnemployed
2                   0.002581339                      0.000000000
3                   0.351691757                      0.004412069
4                   0.425644713                      0.003913733
  employment_changeDecreased employment employment_changeFurloughed
2                          1.764171e-09                0.0005837308
3                          2.352654e-03                0.1642599131
4                          4.654812e-02                0.0035693475
  employment_changeIncreased employment key_workerGovernment defined key worker
2                           0.003158896                              0.04115641
3                           0.615542229                              0.02567268
4                           0.939310952                              0.04082496
  anxiety_disordersAnxiety disorder depressive_disordersDepressive disorder
2                                 0                            0.000000e+00
3                                 0                            5.885981e-11
4                                 0                            8.881784e-16
  eating_disordersEating disorder
2                    0.000000e+00
3                    1.054820e-06
4                    6.354079e-04
  obsessive_compulsive_disordersObsessive compulsive disorder
2                                                2.270939e-11
3                                                3.086507e-04
4                                                9.341355e-04
  psychotic_disordersPsychotic disorder
2                             0.1062995
3                             0.1060261
4                             0.9562884
  personality_disordersPersonality disorder bipolar_disordersBipolar disorder
2                              1.767475e-13                       0.001291487
3                              1.328290e-02                       0.037111202
4                              2.775191e-04                       0.028033458

show p for two tailed z test for ethnicity as is model

p_as.is.gad 
  (Intercept)     age16-18     age19-25     age36-45     age46-55     age56-65
2           0 1.532459e-10 2.678469e-09 4.293519e-05 2.786464e-06 2.220446e-16
3           0 5.674237e-07 2.072432e-04 7.090886e-03 1.092044e-09 3.021583e-12
4           0 6.298919e-01 3.333184e-03 2.215439e-01 9.678825e-01 2.628730e-05
      age66-70   age71-75    age76+ genderFemale
2 4.513545e-06 0.02751186 0.1631106 4.811381e-04
3 6.257205e-04 0.01374216 0.2619255 4.637595e-08
4 1.644757e-02 0.42563666 0.8621599 1.115383e-04
  genderNon-binary/Prefer to self define ethnicityMixed or multiple ethnicity
2                            0.001070268                            0.7716402
3                            0.367655376                            0.5218096
4                            0.036936476                            0.7202514
  ethnicityAsian ethnicityBlack ethnicityArab ethnicityOther
2      0.9990833      0.7594967     0.3532929     0.11238805
3      0.5930515      0.4832388     0.9240917     0.01288128
4      0.0241364      0.1256650     0.9229559     0.31736037
  employment_prior_covidRetired employment_prior_covidStudent
2                   0.003525701                   0.002441345
3                   0.595042224                   0.339955648
4                   0.210032339                   0.418381579
  employment_prior_covidUnemployed employment_changeDecreased employment
2                      0.000000000                          2.119847e-09
3                      0.004749336                          2.444696e-03
4                      0.003923167                          4.532384e-02
  employment_changeFurloughed employment_changeIncreased employment
2                0.0005822465                            0.00326245
3                0.1696790547                            0.64066585
4                0.0031593431                            0.93788200
  key_workerGovernment defined key worker anxiety_disordersAnxiety disorder
2                              0.03974843                                 0
3                              0.02556546                                 0
4                              0.04236972                                 0
  depressive_disordersDepressive disorder eating_disordersEating disorder
2                            0.000000e+00                    0.000000e+00
3                            7.051693e-11                    1.274027e-06
4                            6.661338e-16                    6.861704e-04
  obsessive_compulsive_disordersObsessive compulsive disorder
2                                                2.051292e-11
3                                                2.882067e-04
4                                                8.755800e-04
  psychotic_disordersPsychotic disorder
2                             0.1104281
3                             0.1093510
4                             0.9461185
  personality_disordersPersonality disorder bipolar_disordersBipolar disorder
2                              1.858513e-13                       0.001322596
3                              1.316989e-02                       0.039650249
4                              2.773318e-04                       0.027485825

Relative risk

Relative risk is the ratio of the probability of choosing one outcome category over the probability of choosing the baseline category.

Relative risk of 1 means there is no difference between outcome and baseline. Greater than 1 indicates outcome is more likely than baseline. Less than one indicates it is less likely than baseline.

RR 1.4 would mean there is a 1.4% chance of outcome occurring relative to baseline etc

if something is below one (e.g. 0.6 x the risk) you can also interpret this as 40% less probability of outcomes (1-0.6 = 0.4). So, if women are 0.6 x as likely to belong to a class as men, that means that there is a 40% lower probability of being in that category for women.

exponentiate the model coefficient to obtain relative risk for all predictors. Can see equations here

clustered ethnicity
relative_risk_clustered.gad <- data.frame(exp(coef(gad.model_clustered))) %>%
  rownames_to_column() %>%
  rename(class = "rowname",
    intercept = "X.Intercept.") %>%
  pivot_longer(!class) %>%
  rename(predictor = name,
         RelativeRisk = value) 


relative_risk_clustered.gad$predictor <- gsub("76.","76+",relative_risk_clustered.gad$predictor)


relative_risk_clustered.gad 
ethnicity
relative_risk_as.is.gad <-data.frame(exp(coef(gad.model))) %>%
  rownames_to_column() %>%
  rename(class = "rowname",
    intercept = "X.Intercept.") %>%
  pivot_longer(!class) %>%
  rename(predictor = name,
         RelativeRisk = value)


relative_risk_as.is.gad$predictor <- gsub("76.","76+",relative_risk_as.is.gad$predictor)

head(relative_risk_as.is.gad)

confidence intervals (rr scale)

The confint package has a specific method that is run for multinom objects (see methods(confint)) Obtain the exponent of upper and lower to get relative risk version

First create the CI (converted to relative risk) and reformat it so it is ready to merge with relative risk data for plotting

clustered ethnicity

CI_clustered.gad <- data.frame(exp(confint(gad.model_clustered,level=0.95))) %>%
  rownames_to_column()

names(CI_clustered.gad) <- c("predictor",
                              "2_lower","2_upper",
                              "3_lower","3_upper",
                              "4_lower","4_upper"
                              )

CI_clustered.gad <- CI_clustered.gad %>%
  pivot_longer(!predictor) %>%
  separate(name, c("class","CI")) %>%
  pivot_wider(names_from = CI,
              values_from= value) 


CI_clustered.gad$predictor <- gsub("-",".",CI_clustered.gad$predictor)
CI_clustered.gad$predictor <- gsub(" ",".",CI_clustered.gad$predictor)
CI_clustered.gad$predictor <- gsub("/",".",CI_clustered.gad$predictor)
CI_clustered.gad$predictor <- gsub("\\(Intercept\\)","intercept",CI_clustered.gad$predictor)


  

head(CI_clustered.gad)

ethnicity as is

CI_as.is.gad<- data.frame(exp(confint(gad.model,level=0.95))) %>%
  rownames_to_column()

names(CI_as.is.gad) <-  c("predictor",
                              "2_lower","2_upper",
                              "3_lower","3_upper",
                              "4_lower","4_upper"
                              )

CI_as.is.gad <- CI_as.is.gad %>%
  pivot_longer(!predictor) %>%
  separate(name, c("class","CI")) %>%
  pivot_wider(names_from = CI,
              values_from= value)

CI_as.is.gad$predictor <- gsub("-",".",CI_as.is.gad$predictor)
CI_as.is.gad$predictor <- gsub(" ",".",CI_as.is.gad$predictor)
CI_as.is.gad$predictor <- gsub("/",".",CI_as.is.gad$predictor)
CI_as.is.gad$predictor <- gsub("\\(Intercept\\)","intercept",CI_as.is.gad$predictor)

head(CI_as.is.gad)

forest plots of relative risk

GAD

prepare and format plot data

plot_gad_clustered_dat$category <-plot_gad_clustered_dat$predictor

plot_gad_clustered_dat <-
  plot_gad_clustered_dat %>%
  mutate(category = case_when(grepl("age",predictor) ~ "Age",
                              grepl("gender",predictor) ~ "Gender",
                              grepl("ethnicity_clustered",predictor) ~ "Ethnicity",
                              grepl("employment_prior_covid",predictor) ~ "Employment before covid",
                              grepl("employment_change",predictor) ~ "Change in employment due to covid",
                              grepl("key_worker",predictor) ~ "Key worker status",
                              grepl("anxiety_disorders",predictor) ~ "Pre-existing mental health diagnosis",
                              grepl("depressive_disorders",predictor) ~ "Pre-existing mental health diagnosis",
                              grepl("eating_disorders",predictor) ~ "Pre-existing mental health diagnosis",
                              grepl("obsessive_compulsive_disorders",predictor) ~ "Pre-existing mental health diagnosis",
                              grepl("psychotic_disorders",predictor) ~ "Pre-existing mental health diagnosis",
                              grepl("personality_disorders",predictor) ~ "Pre-existing mental health diagnosis",
                              grepl("bipolar_disorders",predictor) ~ "Pre-existing mental health diagnosis",
                          
                              grepl("intercept",predictor) ~ "Intercept"
                              ),
         
         predictor = case_when(grepl("age",predictor) ~ gsub("age","",predictor),
                               grepl("gender",predictor) ~ gsub("gender","",predictor),
                               grepl("ethnicity_clustered",predictor) ~ gsub("ethnicity_clustered","",predictor),
                               grepl("employment_prior_covid",predictor) ~ gsub("employment_prior_covid","",predictor),
                               grepl("employment_change",predictor) ~ gsub("employment_change","",predictor),
                               grepl("key_worker",predictor) ~ gsub("key_worker","",predictor),
                               grepl("anxiety_disorders",predictor) ~ gsub("anxiety_disorders","",predictor),
                               grepl("depressive_disorders",predictor) ~ gsub("depressive_disorders","",predictor),
                               grepl("eating_disorders",predictor) ~ gsub("eating_disorders","",predictor),
                               grepl("obsessive_compulsive_disorders",predictor) ~ gsub("obsessive_compulsive_disorders","",predictor),
                               grepl("psychotic_disorders",predictor) ~ gsub("psychotic_disorders","",predictor),
                               grepl("personality_disorders",predictor) ~ gsub("personality_disorders","",predictor),
                               grepl("bipolar_disorders",predictor) ~ gsub("bipolar_disorders","",predictor),
                               grepl("intercept",predictor) ~ "intercept")
                               ) 
         

plot_gad_as.is_dat$category <-plot_gad_as.is_dat$predictor

plot_gad_as.is_dat <-
plot_gad_as.is_dat %>%
  mutate(category = case_when(grepl("age",predictor) ~ "Age",
                              grepl("gender",predictor) ~ "Gender",
                              grepl("ethnicity",predictor) ~ "Ethnicity",
                              grepl("employment_prior_covid",predictor) ~ "Employment before covid",
                              grepl("employment_change",predictor) ~ "Change in employment due to covid",
                              grepl("key_worker",predictor) ~ "Key worker status",
                              grepl("anxiety_disorders",predictor) ~ "Pre-existing mental health diagnosis",
                              grepl("depressive_disorders",predictor) ~ "Pre-existing mental health diagnosis",
                              grepl("eating_disorders",predictor) ~ "Pre-existing mental health diagnosis",
                              grepl("obsessive_compulsive_disorders",predictor) ~ "Pre-existing mental health diagnosis",
                              grepl("psychotic_disorders",predictor) ~ "Pre-existing mental health diagnosis",
                              grepl("personality_disorders",predictor) ~ "Pre-existing mental health diagnosis",
                              grepl("bipolar_disorders",predictor) ~ "Pre-existing mental health diagnosis",
                               grepl("intercept",predictor) ~ "Intercept"
                              ),
         
         predictor = case_when(grepl("age",predictor) ~ gsub("age","",predictor),
                               grepl("gender",predictor) ~ gsub("gender","",predictor),
                               grepl("ethnicity",predictor) ~ gsub("ethnicity_clustered","",predictor),
                               grepl("employment_prior_covid",predictor) ~ gsub("employment_prior_covid","",predictor),
                               grepl("employment_change",predictor) ~ gsub("employment_change","",predictor),
                               grepl("key_worker",predictor) ~ gsub("key_worker","",predictor),
                               grepl("anxiety_disorders",predictor) ~ gsub("anxiety_disorders","",predictor),
                               grepl("depressive_disorders",predictor) ~ gsub("depressive_disorders","",predictor),
                               grepl("eating_disorders",predictor) ~ gsub("eating_disorders","",predictor),
                               grepl("obsessive_compulsive_disorders",predictor) ~ gsub("obsessive_compulsive_disorders","",predictor),
                               grepl("psychotic_disorders",predictor) ~ gsub("psychotic_disorders","",predictor),
                               grepl("personality_disorders",predictor) ~ gsub("personality_disorders","",predictor),
                               grepl("bipolar_disorders",predictor) ~  gsub("mental_health_controls","",predictor),
                               grepl("intercept",predictor) ~ "intercept")
                               )

sort variables so categories are together

plot_gad_clustered_dat$category <- factor(plot_gad_clustered_dat$category,
                                       levels=c("Intercept","Gender","Ethnicity","Age","Employment before covid","Change in employment due to covid","Key worker status","Pre-existing mental health diagnosis"))

plot_gad_clustered_dat$predictor <- fct_reorder2(plot_gad_clustered_dat$predictor, plot_gad_clustered_dat$class,  plot_gad_clustered_dat$category)




plot_gad_as.is_dat$category <- factor(plot_gad_as.is_dat$category,
                                       levels=c("Intercept","Gender","Ethnicity","Age","Employment before covid","Change in employment due to covid","Key worker status","Pre-existing mental health diagnosis"))

plot_gad_as.is_dat$predictor <- fct_reorder2(plot_gad_as.is_dat$predictor, plot_gad_as.is_dat$class,  plot_gad_as.is_dat$category)

forest plot clustered ethnicity

Note: smaller category error bars extend beyond a reasonable axis. Wil need fixing for publication.

gad_eth_clustered_forest <-relrisk.forest(plot_gad_clustered_dat,
plot_title="Relative risk (95% CI) of class membership: anxiety",
sub_title="ethnicity clustered",
lowlim = 0.001,
uplim = 10)


gad_eth_clustered_forest
Warning: Removed 1 rows containing missing values (geom_point).

forest plot ethnicity as is

Note: smaller category error bars extend beyond a reasonable axis. Wil need fixing for publication.

gad_eth_as.is_forest <- relrisk.forest(plot_gad_as.is_dat,
plot_title="Relative risk (95% CI) of class membership: anxiety",
sub_title="ethnicity as is",
lowlim = 0.001,
uplim = 10)

gad_eth_as.is_forest
Warning: Removed 2 rows containing missing values (geom_point).

#### save plots

ggsave(
  filename=file.path(dirname(dirname(getwd())),"output/ForestPlots/GAD/anxiety_RelativeRisk_forestPlot_ethnicityClustered.png"),
  plot=gad_eth_clustered_forest, 
  width = 18, height = 10, dpi = 150, units = "in")
Warning: Removed 1 rows containing missing values (geom_point).
ggsave(
  filename=file.path(dirname(dirname(getwd())),"output/ForestPlots/GAD/anxiety_RelativeRisk_forestPlot.png"),
  plot=gad_eth_as.is_forest, 
  width = 18, height = 10, dpi = 150, units = "in")
Warning: Removed 2 rows containing missing values (geom_point).

MASQ

prepare and run models

specify reference category for outcome

will use the largest class group as reference

BiggestClass.masq <- masq.AllData %>%
  group_by(MostLikelyClass) %>%
  summarise(n=n()) %>%
  arrange(-n) %>%
  slice(1)%>%
  select(MostLikelyClass) %>%
  as.character(.)

BiggestClass.masq
[1] "3"

Will need to adapt this to the correct number of classes depending on final model

masq.AllData$MostLikelyClass <- factor(masq.AllData$MostLikelyClass,
                                       levels=c("1","2","3","4","5")) %>%
  relevel(masq.AllData$MostLikelyClass, ref = BiggestClass.masq)
Perform regression (ethnicity clustered)
masq.model_clustered <- multinom(formula_clustered,
                            data = masq.AllData)
# weights:  135 (104 variable)
initial  value 22086.316472 
iter  10 value 12682.548064
iter  20 value 11859.394528
iter  30 value 11165.449558
iter  40 value 10972.460273
iter  50 value 10953.018144
iter  60 value 10949.412688
iter  70 value 10948.984909
iter  80 value 10948.782059
iter  90 value 10948.621659
iter 100 value 10948.568540
final  value 10948.568540 
stopped after 100 iterations
Perform regression (ethnicity as is)
masq.model<- multinom(formula,
                            data = masq.AllData)
# weights:  155 (120 variable)
initial  value 22086.316472 
iter  10 value 12681.333588
iter  20 value 11856.224058
iter  30 value 11165.483136
iter  40 value 10967.798201
iter  50 value 10947.912165
iter  60 value 10944.405467
iter  70 value 10943.778566
iter  80 value 10943.480845
iter  90 value 10943.189106
iter 100 value 10943.069793
final  value 10943.069793 
stopped after 100 iterations

View and process results

1: coefficients 2: std errors 3: residual deviances & AIC ##### view summary: Ethnicity clustered

summary(masq.model_clustered)
Call:
multinom(formula = formula_clustered, data = masq.AllData)

Coefficients:
  (Intercept)   age16-18    age19-25   age36-45    age46-55   age56-65
1  -3.3534339 -0.2246741 -0.27010497 -0.0480830 -0.44509283 -0.4155389
2   0.1708474 -0.3079949  0.07977917 -0.1449227 -0.25211813 -0.1401873
4  -3.4052826 -0.7173888  0.11844416 -0.1225591 -0.02773374 -0.2478546
5  -3.7150027  0.3529555 -0.08392221 -0.3424730 -0.33277432 -0.2616829
     age66-70    age71-75     age76+ genderFemale
1 -0.53716100 -0.56984866 -5.9507835    0.1569743
2  0.09093635 -0.03069305  0.5692933   -0.4114874
4 -0.70999920 -0.36977132 -5.8712138    0.5539187
5 -0.86858283 -0.67315625  1.7192703    0.4443630
  genderNon-binary/Prefer to self define
1                            -0.06891821
2                            -0.29824258
4                            -8.79013616
5                             0.53966418
  ethnicity_clusteredMinoritised ethnic community employment_prior_covidRetired
1                                      0.09033541                    -0.6183422
2                                      0.21639694                    -0.4655660
4                                      0.17409120                    -0.6511315
5                                      0.21466303                    -0.7292526
  employment_prior_covidStudent employment_prior_covidUnemployed
1                     0.8283103                       -0.5044174
2                    -0.7928715                       -0.8130652
4                    -6.0582890                       -0.8721267
5                     0.9471081                       -0.3800778
  employment_changeDecreased employment employment_changeFurloughed
1                           -0.02202836                 -0.16141492
2                           -0.14918425                  0.01364417
4                            0.28898275                  0.24859957
5                            0.55499016                  0.40275886
  employment_changeIncreased employment key_workerGovernment defined key worker
1                           -0.28851977                             -0.12223381
2                            0.05326513                              0.06049872
4                            0.01206147                              0.09970771
5                            0.27626513                              0.00561951
  anxiety_disordersAnxiety disorder depressive_disordersDepressive disorder
1                       -0.04550084                               0.1915367
2                       -0.45641233                              -0.9242325
4                       -0.01059467                              -0.5770875
5                       -0.23254922                               0.1814364
  eating_disordersEating disorder
1                      0.02221847
2                     -0.35745658
4                     -0.09673019
5                      0.11384672
  obsessive_compulsive_disordersObsessive compulsive disorder
1                                                 -0.49238729
2                                                 -0.01646093
4                                                 -0.07723958
5                                                 -0.70983002
  psychotic_disordersPsychotic disorder
1                            -1.1022739
2                             0.3006774
4                             1.0650131
5                             0.3958391
  personality_disordersPersonality disorder bipolar_disordersBipolar disorder
1                                -0.1030690                         0.4655088
2                                -0.6311304                        -0.4095350
4                                -0.5666792                        -0.4232159
5                                -0.0612022                         0.1286832

Std. Errors:
  (Intercept)  age16-18   age19-25   age36-45   age46-55   age56-65  age66-70
1  0.22793179 0.6044032 0.26424974 0.17419140 0.18202038 0.20970425 0.5344713
2  0.07231369 0.2523324 0.09883712 0.06933067 0.06541695 0.06937663 0.1287773
4  0.21660033 0.7301498 0.23101953 0.17745715 0.16409047 0.19202904 0.5300329
5  0.23057104 0.4156004 0.23094755 0.18222530 0.17166663 0.19391529 0.6023694
   age71-75     age76+ genderFemale genderNon-binary/Prefer to self define
1 0.7356971 24.6320715   0.16184059                              0.6101184
2 0.1702867  0.3317533   0.04685269                              0.2310442
4 0.6095128 24.4099826   0.15734105                             53.6686568
5 0.7317874  0.6572730   0.16710814                              0.4910483
  ethnicity_clusteredMinoritised ethnic community employment_prior_covidRetired
1                                       0.3145273                     1.0338619
2                                       0.1108999                     0.2397209
4                                       0.2899838                     1.0302742
5                                       0.2917672                     1.0420769
  employment_prior_covidStudent employment_prior_covidUnemployed
1                     0.7549094                        0.2643221
2                     0.7501142                        0.1118174
4                    25.5920219                        0.3273417
5                     0.7570266                        0.2643364
  employment_changeDecreased employment employment_changeFurloughed
1                            0.17899055                   0.1882251
2                            0.06403244                   0.0611963
4                            0.15679799                   0.1581659
5                            0.15693288                   0.1657373
  employment_changeIncreased employment key_workerGovernment defined key worker
1                            0.30370305                              0.13615317
2                            0.09680611                              0.04581351
4                            0.26006389                              0.12367543
5                            0.25011468                              0.13064332
  anxiety_disordersAnxiety disorder depressive_disordersDepressive disorder
1                        0.16027616                               0.1680468
2                        0.05564113                               0.0535519
4                        0.14619694                               0.1433618
5                        0.15128673                               0.1574589
  eating_disordersEating disorder
1                       0.2164463
2                       0.1009802
4                       0.2185169
5                       0.1992213
  obsessive_compulsive_disordersObsessive compulsive disorder
1                                                  0.26327149
2                                                  0.09771346
4                                                  0.23141461
5                                                  0.26974454
  psychotic_disordersPsychotic disorder
1                             0.7281099
2                             0.2010672
4                             0.3595626
5                             0.3902274
  personality_disordersPersonality disorder bipolar_disordersBipolar disorder
1                                 0.3059850                         0.2834880
2                                 0.1721049                         0.1649833
4                                 0.3804970                         0.3885685
5                                 0.2897254                         0.2971946

Residual Deviance: 21897.14 
AIC: 22105.14 
view summary: Ethnicity as is
summary(masq.model)
Call:
multinom(formula = formula, data = masq.AllData)

Coefficients:
  (Intercept)   age16-18    age19-25    age36-45    age46-55   age56-65
1  -3.3532860 -0.2186131 -0.26837948 -0.04825197 -0.44422205 -0.4145852
2   0.1693906 -0.3054205  0.08206832 -0.14505541 -0.25111543 -0.1394304
4  -3.4157541 -0.7300027  0.11756126 -0.12474580 -0.02557965 -0.2450459
5  -3.7277774  0.3376130 -0.08507963 -0.34766775 -0.32893789 -0.2540240
     age66-70    age71-75     age76+ genderFemale
1 -0.53481212 -0.57099237 -6.4180629    0.1590440
2  0.09195194 -0.02995838  0.5719441   -0.4113503
4 -0.70591381 -0.36205067 -6.3226809    0.5548626
5 -0.88268331 -0.66980771  1.7208891    0.4498306
  genderNon-binary/Prefer to self define ethnicityMixed or multiple ethnicity
1                            -0.05986941                           0.01169794
2                            -0.30056502                           0.25431192
4                            -9.87255482                           0.29701758
5                             0.55041955                          -0.02701734
  ethnicityAsian ethnicityBlack ethnicityArab ethnicityOther
1      0.1473725    0.001870045     -2.859050     0.30483369
2      0.2648712    0.149374753     -6.025050     0.05759361
4      0.5142488   -0.237826749     -3.934510    -8.14998843
5      0.8821718   -7.602556799     -3.199124     0.21387162
  employment_prior_covidRetired employment_prior_covidStudent
1                    -0.6161812                     0.8267753
2                    -0.4694438                    -0.7942560
4                    -0.6532973                    -6.5678224
5                    -0.7278842                     0.9554902
  employment_prior_covidUnemployed employment_changeDecreased employment
1                       -0.5077002                           -0.02490025
2                       -0.8139698                           -0.14776472
4                       -0.8637742                            0.29550063
5                       -0.3812142                            0.55818983
  employment_changeFurloughed employment_changeIncreased employment
1                 -0.16366287                           -0.29223232
2                  0.01388549                            0.05460684
4                  0.25515122                            0.02355797
5                  0.41135488                            0.28909574
  key_workerGovernment defined key worker anxiety_disordersAnxiety disorder
1                            -0.124341183                      -0.044272612
2                             0.060472006                      -0.454657977
4                             0.101148767                      -0.002542874
5                             0.001245017                      -0.223412441
  depressive_disordersDepressive disorder eating_disordersEating disorder
1                               0.1902082                      0.02170131
2                              -0.9250352                     -0.35626486
4                              -0.5760062                     -0.08783083
5                               0.1813459                      0.11619524
  obsessive_compulsive_disordersObsessive compulsive disorder
1                                                 -0.49361826
2                                                 -0.01812124
4                                                 -0.08471810
5                                                 -0.71054890
  psychotic_disordersPsychotic disorder
1                            -1.1089487
2                             0.2979795
4                             1.0585645
5                             0.3955389
  personality_disordersPersonality disorder bipolar_disordersBipolar disorder
1                               -0.10297512                         0.4668157
2                               -0.63152437                        -0.4062939
4                               -0.57089686                        -0.4158154
5                               -0.06450619                         0.1346350

Std. Errors:
  (Intercept)  age16-18   age19-25   age36-45   age46-55   age56-65  age66-70
1  0.22811500 0.6032110 0.26430404 0.17427128 0.18209485 0.20977210 0.5336435
2  0.07236235 0.2520531 0.09886115 0.06934149 0.06542848 0.06939268 0.1287891
4  0.21686631 0.7315828 0.23115207 0.17747905 0.16413489 0.19207318 0.5301621
5  0.23077704 0.4167793 0.23100012 0.18229206 0.17175203 0.19401628 0.6060532
   age71-75     age76+ genderFemale genderNon-binary/Prefer to self define
1 0.7363648 31.1188540    0.1619219                              0.6082385
2 0.1702998  0.3316927    0.0468769                              0.2313392
4 0.6085057 30.7023991    0.1574121                             92.8201235
5 0.7327581  0.6590994    0.1671625                              0.4921610
  ethnicityMixed or multiple ethnicity ethnicityAsian ethnicityBlack
1                            0.4604446      0.5948947      1.0195040
2                            0.1608439      0.1915771      0.3384476
4                            0.3929381      0.4669255      1.0197285
5                            0.4604956      0.4301512     41.8687034
  ethnicityArab ethnicityOther employment_prior_covidRetired
1      17.15678      0.7271190                     1.0324222
2      25.91137      0.3465126                     0.2398563
4      20.93417     49.9348290                     1.0302115
5      18.21465      0.7273932                     1.0453848
  employment_prior_covidStudent employment_prior_covidUnemployed
1                     0.7558608                        0.2643839
2                     0.7501030                        0.1118624
4                    32.9537476                        0.3271549
5                     0.7564152                        0.2641465
  employment_changeDecreased employment employment_changeFurloughed
1                            0.17912629                  0.18834924
2                            0.06405068                  0.06121457
4                            0.15679671                  0.15822370
5                            0.15705322                  0.16594238
  employment_changeIncreased employment key_workerGovernment defined key worker
1                            0.30415164                              0.13619318
2                            0.09683825                              0.04581928
4                            0.25971806                              0.12366495
5                            0.25016794                              0.13069846
  anxiety_disordersAnxiety disorder depressive_disordersDepressive disorder
1                         0.1604032                              0.16814462
2                         0.0556844                              0.05358782
4                         0.1462790                              0.14340352
5                         0.1515374                              0.15768240
  eating_disordersEating disorder
1                       0.2164125
2                       0.1009938
4                       0.2186787
5                       0.1993455
  obsessive_compulsive_disordersObsessive compulsive disorder
1                                                  0.26337714
2                                                  0.09774627
4                                                  0.23177739
5                                                  0.26986582
  psychotic_disordersPsychotic disorder
1                             0.7303503
2                             0.2011655
4                             0.3593871
5                             0.3904275
  personality_disordersPersonality disorder bipolar_disordersBipolar disorder
1                                 0.3060431                         0.2834739
2                                 0.1721302                         0.1649351
4                                 0.3802417                         0.3883717
5                                 0.2896535                         0.2976452

Residual Deviance: 21886.14 
AIC: 22126.14 
create z and p values

make Z

z_clustered.masq <- summary(masq.model_clustered)$coefficients/summary(masq.model_clustered)$standard.errors
  
z_as.is.masq <- summary(masq.model)$coefficients/summary(masq.model)$standard.errors

print z clusted ethnicity model

z_clustered.masq
  (Intercept)   age16-18   age19-25   age36-45   age46-55  age56-65   age66-70
1  -14.712445 -0.3717288 -1.0221580 -0.2760355 -2.4452911 -1.981547 -1.0050325
2    2.362587 -1.2205919  0.8071782 -2.0903119 -3.8540186 -2.020670  0.7061521
4  -15.721502 -0.9825227  0.5127019 -0.6906407 -0.1690149 -1.290714 -1.3395380
5  -16.112182  0.8492666 -0.3633821 -1.8793933 -1.9384915 -1.349470 -1.4419439
    age71-75     age76+ genderFemale genderNon-binary/Prefer to self define
1 -0.7745697 -0.2415868    0.9699317                             -0.1129587
2 -0.1802433  1.7160142   -8.7825768                             -1.2908465
4 -0.6066670 -0.2405251    3.5204971                             -0.1637853
5 -0.9198796  2.6157630    2.6591341                              1.0990044
  ethnicity_clusteredMinoritised ethnic community employment_prior_covidRetired
1                                       0.2872101                    -0.5980897
2                                       1.9512819                    -1.9421174
4                                       0.6003481                    -0.6319982
5                                       0.7357340                    -0.6998069
  employment_prior_covidStudent employment_prior_covidUnemployed
1                     1.0972313                        -1.908343
2                    -1.0570011                        -7.271362
4                    -0.2367257                        -2.664270
5                     1.2510895                        -1.437856
  employment_changeDecreased employment employment_changeFurloughed
1                             -0.123070                  -0.8575631
2                             -2.329823                   0.2229575
4                              1.843026                   1.5717644
5                              3.536481                   2.4301045
  employment_changeIncreased employment key_workerGovernment defined key worker
1                           -0.95000616                             -0.89776690
2                            0.55022491                              1.32054326
4                            0.04637887                              0.80620465
5                            1.10455385                              0.04301414
  anxiety_disordersAnxiety disorder depressive_disordersDepressive disorder
1                       -0.28389025                                1.139782
2                       -8.20278652                              -17.258634
4                       -0.07246849                               -4.025391
5                       -1.53714224                                1.152278
  eating_disordersEating disorder
1                       0.1026512
2                      -3.5398664
4                      -0.4426668
5                       0.5714587
  obsessive_compulsive_disordersObsessive compulsive disorder
1                                                  -1.8702644
2                                                  -0.1684613
4                                                  -0.3337714
5                                                  -2.6314898
  psychotic_disordersPsychotic disorder
1                             -1.513884
2                              1.495407
4                              2.961968
5                              1.014381
  personality_disordersPersonality disorder bipolar_disordersBipolar disorder
1                                -0.3368433                         1.6420758
2                                -3.6671259                        -2.4822820
4                                -1.4893134                        -1.0891670
5                                -0.2112421                         0.4329931

print z ethnicity as is

z_as.is.masq
  (Intercept)   age16-18   age19-25   age36-45   age46-55  age56-65   age66-70
1  -14.699980 -0.3624156 -1.0154195 -0.2768785 -2.4395091 -1.976360 -1.0021899
2    2.340867 -1.2117310  0.8301373 -2.0918993 -3.8380139 -2.009296  0.7139727
4  -15.750506 -0.9978401  0.5085884 -0.7028762 -0.1558453 -1.275795 -1.3315055
5  -16.153155  0.8100522 -0.3683099 -1.9072018 -1.9151906 -1.309292 -1.4564451
    age71-75     age76+ genderFemale genderNon-binary/Prefer to self define
1 -0.7754205 -0.2062436    0.9822266                            -0.09843081
2 -0.1759155  1.7243190   -8.7751168                            -1.29923928
4 -0.5949832 -0.2059344    3.5249038                            -0.10636222
5 -0.9140912  2.6109707    2.6909784                             1.11837290
  ethnicityMixed or multiple ethnicity ethnicityAsian ethnicityBlack
1                           0.02540575      0.2477286    0.001834269
2                           1.58111029      1.3825831    0.441352696
4                           0.75588898      1.1013508   -0.233225567
5                          -0.05867014      2.0508415   -0.181580899
  ethnicityArab ethnicityOther employment_prior_covidRetired
1    -0.1666426      0.4192349                    -0.5968306
2    -0.2325253      0.1662093                    -1.9571874
4    -0.1879468     -0.1632125                    -0.6341390
5    -0.1756347      0.2940247                    -0.6962835
  employment_prior_covidStudent employment_prior_covidUnemployed
1                     1.0938195                        -1.920315
2                    -1.0588626                        -7.276525
4                    -0.1993043                        -2.640260
5                     1.2631821                        -1.443192
  employment_changeDecreased employment employment_changeFurloughed
1                            -0.1390095                  -0.8689330
2                            -2.3069968                   0.2268331
4                             1.8846099                   1.6125979
5                             3.5541444                   2.4789019
  employment_changeIncreased employment key_workerGovernment defined key worker
1                           -0.96081126                             -0.91297661
2                            0.56389739                              1.31979383
4                            0.09070594                              0.81792592
5                            1.15560668                              0.00952587
  anxiety_disordersAnxiety disorder depressive_disordersDepressive disorder
1                       -0.27600821                                1.131218
2                       -8.16490740                              -17.262043
4                       -0.01738372                               -4.016681
5                       -1.47430578                                1.150071
  eating_disordersEating disorder
1                       0.1002775
2                      -3.5275910
4                      -0.4016433
5                       0.5828837
  obsessive_compulsive_disordersObsessive compulsive disorder
1                                                  -1.8741879
2                                                  -0.1853906
4                                                  -0.3655149
5                                                  -2.6329711
  psychotic_disordersPsychotic disorder
1                             -1.518379
2                              1.481265
4                              2.945471
5                              1.013092
  personality_disordersPersonality disorder bipolar_disordersBipolar disorder
1                                -0.3364726                         1.6467679
2                                -3.6688764                        -2.4633563
4                                -1.5014052                        -1.0706633
5                                -0.2227013                         0.4523339

two tailed z test p values

p_clustered.masq <- (1 - pnorm(abs(z_clustered.masq), 0, 1)) * 2
p_as.is.masq <- (1 - pnorm(abs(z_as.is.masq), 0, 1)) * 2

show p for two tailed z test for ethnicity clustered model

p_clustered.masq 
  (Intercept)  age16-18  age19-25   age36-45     age46-55   age56-65  age66-70
1   0.0000000 0.7100948 0.3067061 0.78252084 0.0144735250 0.04752991 0.3148812
2   0.0181479 0.2222406 0.4195638 0.03658980 0.0001161948 0.04331396 0.4800935
4   0.0000000 0.3258424 0.6081598 0.48979138 0.8657849058 0.19680274 0.1803956
5   0.0000000 0.3957330 0.7163195 0.06019081 0.0525632835 0.17718596 0.1493182
   age71-75      age76+ genderFemale genderNon-binary/Prefer to self define
1 0.4385940 0.809100347 0.3320805169                              0.9100633
2 0.8569615 0.086159441 0.0000000000                              0.1967569
4 0.5440719 0.809923197 0.0004307388                              0.8699002
5 0.3576357 0.008902833 0.0078341765                              0.2717662
  ethnicity_clusteredMinoritised ethnic community employment_prior_covidRetired
1                                      0.77395146                    0.54978008
2                                      0.05102352                    0.05212289
4                                      0.54827429                    0.52738803
5                                      0.46189260                    0.48404792
  employment_prior_covidStudent employment_prior_covidUnemployed
1                     0.2725403                     5.634685e-02
2                     0.2905111                     3.559375e-13
4                     0.8128696                     7.715559e-03
5                     0.2109018                     1.504748e-01
  employment_changeDecreased employment employment_changeFurloughed
1                          0.9020516974                  0.39113374
2                          0.0198155077                  0.82356860
4                          0.0653252414                  0.11600520
5                          0.0004054955                  0.01509447
  employment_changeIncreased employment key_workerGovernment defined key worker
1                             0.3421091                               0.3693098
2                             0.5821651                               0.1866537
4                             0.9630083                               0.4201249
5                             0.2693530                               0.9656903
  anxiety_disordersAnxiety disorder depressive_disordersDepressive disorder
1                      7.764945e-01                            2.543771e-01
2                      2.220446e-16                            0.000000e+00
4                      9.422291e-01                            5.688071e-05
5                      1.242585e-01                            2.492069e-01
  eating_disordersEating disorder
1                    0.9182398236
2                    0.0004003297
4                    0.6580067871
5                    0.5676887672
  obsessive_compulsive_disordersObsessive compulsive disorder
1                                                 0.061447109
2                                                 0.866220416
4                                                 0.738552065
5                                                 0.008501141
  psychotic_disordersPsychotic disorder
1                           0.130055305
2                           0.134808239
4                           0.003056793
5                           0.310401087
  personality_disordersPersonality disorder bipolar_disordersBipolar disorder
1                              0.7362350333                        0.10057429
2                              0.0002452921                        0.01305439
4                              0.1364048624                        0.27608025
5                              0.8326983664                        0.66501979

show p for two tailed z test for ethnicity as is model

p_as.is.masq 
  (Intercept)  age16-18  age19-25   age36-45     age46-55   age56-65  age66-70
1  0.00000000 0.7170415 0.3099059 0.78187339 0.0147072314 0.04811399 0.3162519
2  0.01923902 0.2256154 0.4064612 0.03644753 0.0001240335 0.04450577 0.4752440
4  0.00000000 0.3183569 0.6110408 0.48213290 0.8761549532 0.20202808 0.1830227
5  0.00000000 0.4179102 0.7126422 0.05649446 0.0554681961 0.19043538 0.1452696
   age71-75      age76+ genderFemale genderNon-binary/Prefer to self define
1 0.4380912 0.836600671 0.3259882136                              0.9215902
2 0.8603603 0.084650277 0.0000000000                              0.1938618
4 0.5518547 0.836842135 0.0004236366                              0.9152950
5 0.3606689 0.009028562 0.0071242804                              0.2634078
  ethnicityMixed or multiple ethnicity ethnicityAsian ethnicityBlack
1                            0.9797313     0.80434436      0.9985365
2                            0.1138528     0.16679272      0.6589577
4                            0.4497158     0.27074399      0.8155863
5                            0.9532148     0.04028238      0.8559116
  ethnicityArab ethnicityOther employment_prior_covidRetired
1     0.8676512      0.6750444                    0.55062046
2     0.8161300      0.8679922                    0.05032544
4     0.8509184      0.8703511                    0.52599011
5     0.8605809      0.7687390                    0.48625130
  employment_prior_covidStudent employment_prior_covidUnemployed
1                     0.2740341                     5.481815e-02
2                     0.2896624                     3.426148e-13
4                     0.8420247                     8.284237e-03
5                     0.2065237                     1.489663e-01
  employment_changeDecreased employment employment_changeFurloughed
1                           0.889442647                  0.38488379
2                           0.021054999                  0.82055350
4                           0.059482524                  0.10683189
5                           0.000379211                  0.01317875
  employment_changeIncreased employment key_workerGovernment defined key worker
1                             0.3366471                               0.3612548
2                             0.5728240                               0.1869039
4                             0.9277263                               0.4133995
5                             0.2478421                               0.9923996
  anxiety_disordersAnxiety disorder depressive_disordersDepressive disorder
1                      7.825418e-01                            2.579633e-01
2                      2.220446e-16                            0.000000e+00
4                      9.861305e-01                            5.902345e-05
5                      1.403993e-01                            2.501147e-01
  eating_disordersEating disorder
1                    0.9201239860
2                    0.0004193595
4                    0.6879465950
5                    0.5599716010
  obsessive_compulsive_disordersObsessive compulsive disorder
1                                                 0.060904525
2                                                 0.852922650
4                                                 0.714727063
5                                                 0.008464156
  psychotic_disordersPsychotic disorder
1                           0.128918841
2                           0.138535850
4                           0.003224632
5                           0.311016405
  personality_disordersPersonality disorder bipolar_disordersBipolar disorder
1                              0.7365145044                        0.09960576
2                              0.0002436188                        0.01376431
4                              0.1332507900                        0.28432086
5                              0.8237680106                        0.65102845

Relative risk

Relative risk is the ratio of the probability of choosing one outcome category over the probability of choosing the baseline category.

Relative risk of 1 means there is no difference between outcome and baseline. Greater than 1 indicates outcome is more likely than baseline. Less than one indicates it is less likely than baseline.

RR 1.4 would mean there is a 1.4% chance of outcome occurring relative to baseline etc

if something is below one (e.g. 0.6 x the risk) you can also interpret this as 40% less probability of outcomes (1-0.6 = 0.4). So, if women are 0.6 x as likely to belong to a class as men, that means that there is a 40% lower probability of being in that category for women.

exponentiate the model coefficient to obtain relative risk for all predictors. Can see equations here

clustered ethnicity
relative_risk_clustered.masq <- data.frame(exp(coef(masq.model_clustered))) %>%
  rownames_to_column() %>%
  rename(class = "rowname",
    intercept = "X.Intercept.") %>%
  pivot_longer(!class) %>%
  rename(predictor = name,
         RelativeRisk = value) 


relative_risk_clustered.masq$predictor <- gsub("76.","76+",relative_risk_clustered.masq$predictor)


relative_risk_clustered.masq 
ethnicity
relative_risk_as.is.masq <-data.frame(exp(coef(masq.model))) %>%
  rownames_to_column() %>%
  rename(class = "rowname",
    intercept = "X.Intercept.") %>%
  pivot_longer(!class) %>%
  rename(predictor = name,
         RelativeRisk = value)


relative_risk_as.is.masq$predictor <- gsub("76.","76+",relative_risk_as.is.masq$predictor)

head(relative_risk_as.is.masq)

confidence intervals (rr scale)

The confint package has a specific method that is run for multinom objects (see methods(confint)) Obtain the exponent of upper and lower to get relative risk version

First create the CI (converted to relative risk) and reformat it so it is ready to merge with relative risk data for plotting

clustered ethnicity

CI_clustered.masq <- data.frame(exp(confint(masq.model_clustered,level=0.95))) %>%
  rownames_to_column()

names(CI_clustered.masq) <- c("predictor",
                              "1_lower","1_upper",
                              "2_lower","2_upper",
                              "4_lower","4_upper",
                              "5_lower","5_upper"
                              )

CI_clustered.masq <- CI_clustered.masq %>%
  pivot_longer(!predictor) %>%
  separate(name, c("class","CI")) %>%
  pivot_wider(names_from = CI,
              values_from= value) 


CI_clustered.masq$predictor <- gsub("-",".",CI_clustered.masq$predictor)
CI_clustered.masq$predictor <- gsub(" ",".",CI_clustered.masq$predictor)
CI_clustered.masq$predictor <- gsub("/",".",CI_clustered.masq$predictor)
CI_clustered.masq$predictor <- gsub("\\(Intercept\\)","intercept",CI_clustered.masq$predictor)


  

head(CI_clustered.masq)

ethnicity as is

CI_as.is.masq<- data.frame(exp(confint(masq.model,level=0.95))) %>%
  rownames_to_column()

names(CI_as.is.masq) <-  c("predictor",
                              "1_lower","1_upper",
                              "2_lower","2_upper",
                              "4_lower","4_upper",
                              "5_lower","5_upper"
                              )

CI_as.is.masq <- CI_as.is.masq %>%
  pivot_longer(!predictor) %>%
  separate(name, c("class","CI")) %>%
  pivot_wider(names_from = CI,
              values_from= value)

CI_as.is.masq$predictor <- gsub("-",".",CI_as.is.masq$predictor)
CI_as.is.masq$predictor <- gsub(" ",".",CI_as.is.masq$predictor)
CI_as.is.masq$predictor <- gsub("/",".",CI_as.is.masq$predictor)
CI_as.is.masq$predictor <- gsub("\\(Intercept\\)","intercept",CI_as.is.masq$predictor)

head(CI_as.is.masq)

forest plots of relative risk

MASQ

prepare and format plot data

plot_masq_clustered_dat$category <-plot_masq_clustered_dat$predictor

plot_masq_clustered_dat <-
  plot_masq_clustered_dat %>%
  mutate(category = case_when(grepl("age",predictor) ~ "Age",
                              grepl("gender",predictor) ~ "Gender",
                              grepl("ethnicity_clustered",predictor) ~ "Ethnicity",
                              grepl("employment_prior_covid",predictor) ~ "Employment before covid",
                              grepl("employment_change",predictor) ~ "Change in employment due to covid",
                              grepl("key_worker",predictor) ~ "Key worker status",
                              grepl("anxiety_disorders",predictor) ~ "Pre-existing mental health diagnosis",
                              grepl("depressive_disorders",predictor) ~ "Pre-existing mental health diagnosis",
                              grepl("eating_disorders",predictor) ~ "Pre-existing mental health diagnosis",
                              grepl("obsessive_compulsive_disorders",predictor) ~ "Pre-existing mental health diagnosis",
                              grepl("psychotic_disorders",predictor) ~ "Pre-existing mental health diagnosis",
                              grepl("personality_disorders",predictor) ~ "Pre-existing mental health diagnosis",
                              grepl("bipolar_disorders",predictor) ~ "Pre-existing mental health diagnosis",
                              grepl("mental_health_controls",predictor) ~ "Pre-existing mental health diagnosis",
                              grepl("intercept",predictor) ~ "Intercept"
                              ),
         
         predictor = case_when(grepl("age",predictor) ~ gsub("age","",predictor),
                               grepl("gender",predictor) ~ gsub("gender","",predictor),
                               grepl("ethnicity_clustered",predictor) ~ gsub("ethnicity_clustered","",predictor),
                               grepl("employment_prior_covid",predictor) ~ gsub("employment_prior_covid","",predictor),
                               grepl("employment_change",predictor) ~ gsub("employment_change","",predictor),
                               grepl("key_worker",predictor) ~ gsub("key_worker","",predictor),
                               grepl("anxiety_disorders",predictor) ~ gsub("anxiety_disorders","",predictor),
                               grepl("depressive_disorders",predictor) ~ gsub("depressive_disorders","",predictor),
                               grepl("eating_disorders",predictor) ~ gsub("eating_disorders","",predictor),
                               grepl("obsessive_compulsive_disorders",predictor) ~ gsub("obsessive_compulsive_disorders","",predictor),
                               grepl("psychotic_disorders",predictor) ~ gsub("psychotic_disorders","",predictor),
                               grepl("personality_disorders",predictor) ~ gsub("personality_disorders","",predictor),
                               grepl("bipolar_disorders",predictor) ~ gsub("bipolar_disorders","",predictor),
                               grepl("mental_health_controls",predictor) ~ gsub("mental_health_controls","",predictor),
                               grepl("intercept",predictor) ~ "intercept")
                               ) 
         

plot_masq_as.is_dat$category <-plot_masq_as.is_dat$predictor

plot_masq_as.is_dat <-
plot_masq_as.is_dat %>%
  mutate(category = case_when(grepl("age",predictor) ~ "Age",
                              grepl("gender",predictor) ~ "Gender",
                              grepl("ethnicity",predictor) ~ "Ethnicity",
                              grepl("employment_prior_covid",predictor) ~ "Employment before covid",
                              grepl("employment_change",predictor) ~ "Change in employment due to covid",
                              grepl("key_worker",predictor) ~ "Key worker status",
                              grepl("anxiety_disorders",predictor) ~ "Pre-existing mental health diagnosis",
                              grepl("depressive_disorders",predictor) ~ "Pre-existing mental health diagnosis",
                              grepl("eating_disorders",predictor) ~ "Pre-existing mental health diagnosis",
                              grepl("obsessive_compulsive_disorders",predictor) ~ "Pre-existing mental health diagnosis",
                              grepl("psychotic_disorders",predictor) ~ "Pre-existing mental health diagnosis",
                              grepl("personality_disorders",predictor) ~ "Pre-existing mental health diagnosis",
                              grepl("bipolar_disorders",predictor) ~ "Pre-existing mental health diagnosis",
                              grepl("mental_health_controls",predictor) ~ "Pre-existing mental health diagnosis",
                              grepl("intercept",predictor) ~ "Intercept"
                              ),
         
         predictor = case_when(grepl("age",predictor) ~ gsub("age","",predictor),
                               grepl("gender",predictor) ~ gsub("gender","",predictor),
                               grepl("ethnicity",predictor) ~ gsub("ethnicity_clustered","",predictor),
                               grepl("employment_prior_covid",predictor) ~ gsub("employment_prior_covid","",predictor),
                               grepl("employment_change",predictor) ~ gsub("employment_change","",predictor),
                               grepl("key_worker",predictor) ~ gsub("key_worker","",predictor),
                               grepl("anxiety_disorders",predictor) ~ gsub("anxiety_disorders","",predictor),
                               grepl("depressive_disorders",predictor) ~ gsub("depressive_disorders","",predictor),
                               grepl("eating_disorders",predictor) ~ gsub("eating_disorders","",predictor),
                               grepl("obsessive_compulsive_disorders",predictor) ~ gsub("obsessive_compulsive_disorders","",predictor),
                               grepl("psychotic_disorders",predictor) ~ gsub("psychotic_disorders","",predictor),
                               grepl("personality_disorders",predictor) ~ gsub("personality_disorders","",predictor),
                               grepl("bipolar_disorders",predictor) ~ gsub("bipolar_disorders","",predictor),
                               grepl("mental_health_controls",predictor) ~ gsub("mental_health_controls","",predictor),
                               grepl("intercept",predictor) ~ "intercept")
                               )

sort variables so categories are together

plot_masq_clustered_dat$category <- factor(plot_masq_clustered_dat$category,
                                       levels=c("Intercept","Gender","Ethnicity","Age","Employment before covid","Change in employment due to covid","Key worker status","Pre-existing mental health diagnosis"))

plot_masq_clustered_dat$predictor <- fct_reorder2(plot_masq_clustered_dat$predictor, plot_masq_clustered_dat$class,  plot_masq_clustered_dat$category)




plot_masq_as.is_dat$category <- factor(plot_masq_as.is_dat$category,
                                       levels=c("Intercept","Gender","Ethnicity","Age","Employment before covid","Change in employment due to covid","Key worker status","Pre-existing mental health diagnosis"))

plot_masq_as.is_dat$predictor <- fct_reorder2(plot_masq_as.is_dat$predictor, plot_masq_as.is_dat$class,  plot_masq_as.is_dat$category)

forest plot clustered ethnicity

Note: smaller category error bars extend beyond a reasonable axis. Wil need fixing for publication.

masq_eth_clustered_forest <-relrisk.forest(plot_masq_clustered_dat,
plot_title="Relative risk (95% CI) of class membership: anhedonia",
sub_title="ethnicity clustered",
lowlim = 0.001,
uplim = 10)


masq_eth_clustered_forest
Warning: Removed 1 rows containing missing values (geom_point).

forest plot ethnicity as is

Note: smaller category error bars extend beyond a reasonable axis. Wil need fixing for publication.

masq_eth_as.is_forest <- relrisk.forest(plot_masq_as.is_dat,
plot_title="Relative risk (95% CI) of class membership: anhedonia",
sub_title="ethnicity as is",
lowlim = 0.001,
uplim = 10)

masq_eth_as.is_forest
Warning: Removed 3 rows containing missing values (geom_point).

#### save plots

ggsave(
  filename=file.path(dirname(dirname(getwd())),"output/ForestPlots/MASQ/anhedonia_RelativeRisk_forestPlot_ethnicityClustered.png"),
  plot=masq_eth_clustered_forest, 
  width = 18, height = 10, dpi = 150, units = "in")
Warning: Removed 1 rows containing missing values (geom_point).
ggsave(
  filename=file.path(dirname(dirname(getwd())),"output/ForestPlots/MASQ/anhedonia_RelativeRisk_forestPlot.png"),
  plot=masq_eth_as.is_forest, 
  width = 18, height = 10, dpi = 150, units = "in")
Warning: Removed 3 rows containing missing values (geom_point).